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ABSTRACT 


This  research  considers  a  cantilever  beam  which  can  move 
axially  in  and  out  of  a  rigid  frictionless  hole  and  is  free  to 
vibrate  laterally  outside  the  hole .  Two  Euler  equations 
describing  the  lateral  and  axial  motion  of  the  beam  are 
presented.  A  transformation  of  coordinates  to  eliminate  the 
moving  boundary,  and  spatial  non  dimensionalization  are  used 
to  transform  the  problem  into  a  system  of  two  coupled  non 
linear  partial  differential  equations  with  a  fixed  domain.  A 
finite  element  formulation  provides  a  numerical  solution  to 
the  problem.  Results  for  some  problems  are  presented. 
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INTRODUCTION 


I . 

A  literature  search  of  the  engineering  journals  shows  that 
an  investigation  of  the  transient  behavior  of  a  cantilever 
beam,  free  to  move  axially  in  a  frictionless  hole  at  its 
'fixed'  end,  has  not  been  undertaken  to  this  date.  In  1979, 
Boresi  and  Salinas  prepared  a  report  for  the  Naval  Sea  Systems 
Command,  that  formulates  the  problem  and  proposes  a  solution 
procedure.  The  report  was  the  result  of  an  interest  in  the 
transient  behavior  of  a  gun  barrel  during  recoil  following 
firing.  [Ref.l] 

Hamilton's  principle  was  used  to  generate  the  governing 
partial  differential  equations  for  axial  and  lateral  motion  of 
the  beam  [Ref.l].  As  a  result  of  axial  motion  of  the  beam, 
the  length  of  the  beam  changes  with  time.  Thus  the  'free'  end 
of  the  cantilever  beam  is  a  moving  boundary  point.  If  the 
beam  is  subjected  to  an  axial  force,  then  the  beam  length, 
that  is  the  location  of  the  'free'  end,  is  an  unknown  which  is 
itself  a  solution  of  the  problem.  This  is  a  'conjugate' 
problem,  wherein  the  boundary  condition  is  a  solution  of  a 
problem  which  can  not  be  solved  until  the  boundary  extent  is 
known.  The  analogy  is  of  a  dog  chasing  its  own  tail,  or  the 
'catch  22'  syndrome.  The  dilemma  is  resolved  by  introducing 
a  coordinate  transformation  which  produces  a  classical  two- 
point  (fixed)  boundary  domain.  The  removal  of  the  moving 
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boundary  is  not  without  expense,  as  the  resulting  governing 
partial  differential  equations  are  significantly  more 
complicated.  Thus  the  complication  of  the  boundary  condition 
has  been  'transferred'  into  the  interior  domain  of  the 
problem.  The  two  equations  governing  axial  and  lateral 
motion,  for  beam  length  and  lateral  motion,  are  both  coupled 
and  nonlinear  if  the  axial  motion  is  not  prescribed. 

Using  the  finite  element  method  over  the  spatial  domain, 
the  two  partial  differential  equations  in  space  and  time,  are 
reduced  to  a  system  of  ordinary  differential  equations  in  time 
only.  That  is,  the  original  initial- (two-point)  boundary 
value  problem  is  transformed  into  a  system  of  initial-value 
problems  for  the  transient  behavior  at  discrete  points  of  the 
system.  These  nonlinear  O.D.E.'s  are  linearized  using  the 
quasi-linearization  technique  of  Bellman  [Ref. 2],  and  then 
solved  by  using  a  fifth  order  Gear'  method  for  stiff  equations 

This  investigation  adds  further  to  the  formulation  of  the 
problem  by  i_he  introduction  of  non  dimensional  variables. 
Additionally,  the  work  also  provides  mathematical  development 
and  details  required  for  the  numerical  solution  of  the 
problem.  Restrictions  and  a  generalization  of  the  problem  are 
also  discussed. 

The  scope  of  the  problem  suggests  a  cautious  two-stage 
investigation.  In  the  first  stage,  the  axial  motion  as  a 
function  of  time  is  prescribed.  The  result  is  the  elimination 
of  the  need  to  solve  the  equation  for  axial  motion.  However, 
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the  equation  can  be  used  to  solve  for  the  axial  force 
directly.  Moreover,  the  remaining  governing  equation  for 
lateral  transient  behavior  is  linear  since  the  'length'  term 
in  the  equation  is  known.  It  is  felt  that  the  first  stage 
investigation,  which  is  the  body  of  this  thesis  activity, 
would  provide  useful  insight  into  the  nature  of  the  problem 
prior  to  undertaking  the  second  stage  investigation.  In  the 
second  stage  investigation,  instead  of  prescribing  the  axial 
motion,  the  axial  force  at  the  sliding  end  is  prescribed.  As 
a  result,  the  equation  for  the  transient  axial  response  needs 
to  be  solved  in  conjunction  with  the  equation  for  transient 
lateral  response,  since  now  the  length  of  the  beam  is  also 
unknown.  The  second  stage  problem  is  formulated  but  not 
solved  here. 
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II.  PROBLEM  FORMULATION 


Consider  the  transient  behavior  of  a  cantilever  beam 
fitted  snugly  into  a  frictionless  hole  as  shown  in  Figure 
(2.1)  .  The  beam  is  free  to  move  axially  and  laterally  when  an 
axial  force  F(t)  is  applied,  or  when  otherwise  an  axial 
displacement  is  imposed.  The  beam's  motion  can  be  d^'ned 
completely  by  its  axial  motion  u(t)  as  a  function  of  time,  and 
its  lateral  motion  U(x,t)  as  a  function  of  both  time  and 
position  along  the  x  axis.  Because  of  inertia,  under  certain 
conditions,  such  as  when  the  axial  force  F(t)  is  a  large 
magnitude  impulse,  the  axial  movement  of  the  beam  may  tend  to 
bend  the  beam  by  beam-column  action  or  compress  the  beam 
axially  by  beam-bar  action.  These  axial  deformation  effects 
are  not  considered  here,  that  is  u'  =0.  Therefore,  it  is 
assumed  that  all  points  along  the  x-axis  of  the  beam 
experience  the  same  axial  motion.  Thus,  the  instantaneous 
length  of  the  beam,  L(t),  serves  to  describe  the  axial  motion 
of  all  points  of  the  beam. 

As  the  beam  moves  axially,  the  length  of  the  beam  outside 
the  hole  at  any  time  t  is  defined  as  L(t).  Because  L(t)  is 
changing  with  time  the  extent  of  the  domain  of  the  problem 
changes  with  time.  It  is  this  changing  domain  that  results  in 
the  coupling  of  the  equations  which  describe  the  lateral  and 
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axial  motion  of  the  beam.  The  changing  domain  is  the  essence 
of  the  problem  and  will  be  discussed  at  length  in  the 
development  that  follows.  This  investigation  will  be 
restricted  to  long  slender  beams,  which  in  this  case  will  be 
beams  for  which  the  length  is  equal  to  or  greater  than  ten 
times  either  of  the  cross-sectional  dimensions.  With  this 
restriction  imposed,  the  Timoshenko  Beam  shear  effects  and 
rotary  inertia,  are  neglected  [Ref.  3] .  However,  as  the  beam 
length  becomes  shorter  these  effects  become  larger  and  loss  of 
accuracy  in  the  solution  is  expected. 


Figure  2.1  Cantilever  Beam  Moving 
Axially  in  a  Frictionless  Hole 
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A.  THE  EULER  EQUATION  OF  LATERAL  MOTION 

Imposing  equilibrium  in  the  lateral  direction  and  using 
small  displacement  theory  results  in  the  Euler  Equation  for 
the  lateral  motion  of  a  beam, 

El  v^fx,  t)  +  p  vtt  (x,  t)  =  p{x,  t) 

t  >  0  a) 

0  <  x  <  L(t) 


where  the  subscripts  t  and  x  denote  partial  differentiation 
with  respect  to  time  and  position,  respectively  and; 

u(x,t)  =  the  lateral  displacement  as  a  function  of  x  and 
t . 

E  =  Young's  modulus  of  elasticity  of  the  beam. 

I  =  moment  of  inertia  of  the  beam  cross-section, 
p  =  the  mass  of  the  beam  per  unit  length  (constant) . 

P  =  the  internally  applied  load  per  unit  length. 


The  fourth  order  Euler  Equation  has  two  essential  (forced) 
boundary  conditions  on  displacement  and  slope  at  the  'fixed' 
left  end. 


v  (0,  t)  =  0 
vx<0, t) =  0 


(2) 
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and  two  natural  boundary  conditions  on  moment  and  shear  force 
at  the  'free'  right  end, 


El  {L,  t )  =  M 
El  \ xxx  (L,  t)  =  P 


(3) 


where  M  and  P  are  the  applied  moment  and  load,  respectively. 
The  homogeneous  boundary  conditions  (M  =  P  =  0)  are  the 
boundary  conditions  considered  here.  However,  a  verification 
of  the  solution  method  is  presented  where  the  non  homogeneous 
boundary  conditions  are  imposed.  The  term  'fixed'  end  is  used 
in  reference  to  the  boundary  located  at  the  left  end  of  the 
beam's  domain  (See  Figure  2.2),  i.e.,  at  x=0.  As  a  result  of 
the  axial  motion,  the  point  on  the  beam  at  this  left  or 
'fixed'  end  is  changing  with  time. 

The  natural  boundary  conditions  at  the  free  end  (Eqs .  3), 
for  moment  and  shear,  occurs  at  the  right  end  point  of  the 
beam  (i.e.,  at  x  =  L)  for  all  time  t.  It  is  the  fact  that  the 
argument  L  in  Equations  (3)  is  changing  with  time  that  makes 
the  natural  boundary  conditions  troublesome.  These  so  called 
moving  boundary  conditions  (or  changing  domain)  will  be 
discussed  later  at  length. 

The  Euler  Equation  for  lateral  motion  is  also  a  second 
order  differential  equation  in  time.  To  obtain  a  solution, 
two  initial  conditions,  one  on  its  lateral  position  \)(x,0), 
and  one  on  its  velocity  \>t(x,0),  along  the  x  axis  will  be 


7 


needed.  These  initial  conditions  will  depend  on  the  specific 
problem  being  solved. 

B.  THE  EULER  EQUATION  OF  AXIAL  MOTION 

If  F(t)  rather  than  L(t)  is  prescribed,  then  a 
differential  equation  defining  L(t)  is  needed.  Again,  using 
principles  of  equilibrium  for  motion  in  the  axial  direction, 
the  following  Euler  equation  for  axial  behavior  is  obtained, 

£<t)  +  vL(L,  t)  -pvi(L,t)]=i_F(t)  (4) 

*PL0  PL0 


Equation  (4)  is  subject  to  the  initial  conditions, 

H  0)  =  L0 
L(  0)  =  L0 


(5) 


where  L0  is  the  initial  length  of  the  beam  at  time  t  =  0.  The 
dot  and  double  dot  above  L  denote  the  first  and  second 
derivatives  with  respect  to  time  respectively,  that  is  the 
velocity  and  acceleration  of  axial  motion. 

Together  Equations  (1)  and  (4)  along  with  their  respective 
boundary  and  initial  conditions  form  a  coupled  and  nonlinear, 
initial-boundary  value  problem.  When  the  force  F(t)  is  known, 
these  coupled  nonlinear  equations  can  be  solved  using  the 
finite  element  method  with  a  linearization  scheme  to  find 
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\)(x,t)  and  L  (t )  .  When  L(t)  is  specified,  Equation  (4)  yields 
F (t)  directly. 

C.  THE  MOVING  BOUNDARY 

In  the  boundary  conditions  described  in  Equations  (3)  ,  the 
beam  length  L(t)  is  a  function  of  time.  Thus  the  boundary 
conditions  are  conditions  on  a  boundary  which  is  moving. 
Graphically  this  is  shown  in  Figure  (2.2).  The  curved 
boundary  of  the  region  of  integration  presents  a  problem.  The 
desire  is  to  remove  the  argument  of  time  varying  length  from 
the  boundary  condition  at  the  free  end.  In  essence,  we  desire 
to  secure  the  boundary.  Graphically  the  boundary  becomes  a 
straight  line  where  previously  it  was  a  curved  line  (See 
Figure  2.3)  .  This  can  be  achieved  by  using  a  coordinate 
transformation  as  shown  in  the  next  section. 


Figure  2.2  Region  of  Integration  for  Equation  (1) 


Figure  2.3  Region  of  Integration  for  Equation  (25) 

D.  THE  COORDINATE  TRANSFORMATION 

New  variables  %(x,t)  and  T)(t)  are  introduced  as  follows: 

%  -  _JL_ 

Ht)  (6) 

T|  =  t 


It  should  be  noted  that  ^  is  a  non-dimensional  variable  with 
respect  to  the  spatial  domain.  The  lateral  deflection  now 
becomes  a  function  of  these  variables  as  shown  below. 


v  (  x(£,t|)  ,  t  (%,  Tp  )  ,  or 

v{  £,  (X,  t)  ,Tj  (x,  t)  ) 
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Considering  the  relations  defined  in  Equations  (6),  the 
following  partial  derivatives  are  obtained. 


_  ~xL 

It  ' 


L 

in  - 

dt 

in 

Bx 


where  t  -  — 


Bl 

Bt 


»  1 


=  0 


(8) 


1 .  Transformation  of  the  Spatial  Fourth  Derivative  of  v 

Considering  Equation  (7) ,  the  transformation  of  the 
spatial  fourth  derivative  on  lateral  displacement,  to  the 

new  coordinate  %  is  accomplished  through  a  series  of 
differentiations  using  the  chain  rule.  The  first 

differentiation  results  in, 

Bx  Bx  B^  +  Bx  dr| 

a^  ~Bllhc  a^a^  <9) 


Following  the  substitution  of  Equations  (8)  into  Equation  (9)  , 
Equation  (1C)  is  obtained. 


v 


X 


(10) 
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After  another  differentiation  with  the  chain  rule  the  second 
spacial  derivative  is  found. 


d2  v 

!hT2 


dr\ 

dx 


(11) 


Again,  using  Equations  (8),  the  second  derivative  is  equal  to, 


iv 
L2  « 


(12) 


Likewise,  the  third  derivative  is, 


V  -  v 

***  L3V 


(13) 


and  finally  the  fourth  derivative  is, 


iv 

L*  «« 


(14) 


2 .  Transformation  of  the  Time  Second  Derivative  of  i) 

The  transformation  of  the  time  second  derivative  on 
lateral  deflection  (or  acceleration)  ,  \)tt  to  the  new 
coordinates  E,  and  T)  is  performed  in  a  similar  fashion  as  was 
the  transformation  of  it  spacial  derivatives.  Once  again. 
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using  Equation  (7)  and  the  chain  rule,  the  following 
expression  for  the  first  time  derivative  is  obtained. 


v  =  _9v  _  + 

dt  dt,  dt  dr)  dt 


Substituting  the  Equation  (8)  values  of  the  partial 
derivatives  into  Equation  (15)  results  in  the  following 
expression  for  Dt, 


-  -^v  + 
~L  i 


Another  time  derivative  using  the  chain  rule  results  in  the 
following  equation  for  \)tt, 


5  > 

-  Jif-iiv 
at  L  l  « 


-h—v  +  v  ^  ^  + 


Again,  using  Equations  (8) , 


-  -ii  and,  in 

at  l  at 
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along  with  the  product  rule  of  differentiation,  we  obtain 


v 


tt 


+ 


v 


a  it 

lT 


"in 


“1 

f  ,  \ 

-4- 

II 

^  / 

-  E— v  +  v 


(19) 


Recalling  that  the  coordinate  transformation  on  time  stated 
that  t=r|,  it  follows  that, 


L{t)  =  L(T|) 


_  dL 
dt  9r| 


(20) 


Now,  using  the  quotient  rule  of  differentiation,  Equation  (19) 
becomes. 


v 


tt 


tt 


+  V 


_____ 


(21) 


Finally,  after  multiplying  and  collecting  like  terms,  Dtt 
becomes , 


=  E2  —  v  + 

S  L2  « 


'  ^LV 


6 


(22) 
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K.  THE  FINAL  EULER  EQUATIONS 

Using  the  transformed  operators,  the  Euler  Equations  are 
rewritten  in  terms  of  the  new  coordinates,  4  an<*  Tl* 

1 .  The  Transformed  Euler  Equation  for  Lateral  Deflection 
Substituting  the  transformed  operators  from  Equations 
(14)  and  (22) ,  into  the  original  Euler  Equation  for  lateral 
deflection, 

+  pvtt  =  P(x,  t)  (1) 


results  in  the  following  Euler  Equation  transformed  to  the  4 
and  rj  coordinates, 


El 

T?V««  +  P 


V  + 

tt 


24 


3-. 


P(4,T|) 


t  L 
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Multiplying  through  by  the  inverse  of  the  coefficient  of 
gives. 


(  .  N 

z 

2 

V  +  ££l 

4444  -§j 

S2 

t 

v«  *  24 

- 

W 

v  -2£iv  -  +  v 

(25) 


0  <  ^  <  1 
0  <  t] 


and  its  boundary  conditions. 


V  ( 0 , T|)  =  0 

Vw(l,Ttf  =  0 

V5(0,T))  =  0 

VW(1,TU  =  0 

(26) 


The  boundary  conditions  are  now  functions  of  ^  over  the  domain 
0  <  4  <  i  ,  in  lieu  of  x  over  the  domain  0  <  x  <  L(t) .  The 
initial  conditions  on  deflection  and  velocity  will  be 
functions  of  4  as  well. 

2 .  The  Transformed  Euler  Equation  for  Axial  Motion 

The  coordinate  transformation  on  the  Euler  Equation  of 
axial  motion  shown  again  here, 

L{t)  +  _i_  [EIVL  (L,t)  -  p Vt  (L,  t)  1  =  i-F(t)  (4) 

2pL0  1  J  PL0 
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results  in  the  transformed  equation. 


/ 


L  + 


1 


d,T|) 


+  V^(1,T1) 


1 


PA> 


F(T[) 


(28) 


subjected  to  the  initial  conditions  in  Equations  (5) . 

3.  Non  dimensionalization  of  the  Lateral  Deflection,  v 
The  purpose  of  the  coordinate  transformation  just 
completed  was  to  deal  with  the  difficulty  presented  by  the 
moving  boundary  condition  at  the  free  end  of  the  beam.  The 
four  boundary  conditions  of  Equations  (2)  and  (3)  were  also 
transformed  to  the  £  and  TJ  coordinates  as  shown  in  Equations 
(26)  .  One  of  the  great  difficulties  encountered  in  this 
investigation  resulted  from  the  coordinate  transformation 
performed  on  the  boundary  conditions.  After  the  introduction 
of  the  non  dimensional  variable  %,  the  finite  element  method 
(FEM)  of  Chapter  III  was  pursued.  This  included  an  attempt  to 
confirm  the  FEM  program  on  a  couple  of  statics  problems  with 
known  solutions.  The  resulting  FEM  solutions  were  L3  larger 
for  displacements  and  L2  larger  for  slopes.  We  had  simply 
imposed  the  load  (or  moment)  as  one  would  have  if  the  problem 
had  the  dimensional  independent  variable  x,  when  in  fact  x  had 
been  replaced  by  the  non  dimensional  variable  £  =  x/L.  This 
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problem  was  eventually  resolved  by  introducing  the  non 
dimensional  displacement,  X>* ,  defined  as, 


v* 


(29) 


Its  derivatives, 


d\‘  _  1 

dv  L 


(30) 


are  used  such  that. 


V*  = 


3v 


—  (LV  )  =  L\l 
=  “  <£VS  >  -  £v« 


(31) 


and  in  the  same  fashion, 


(32) 


and, 


(33) 
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After  making  the  substitutions  into  the  equations  of 
lateral  and  axial  motion.  Equations  (25)  and  (28) 
respectively,  the  spatially  non  dimensional  Euler  equations 
are  obtained. 

a.  Final  Euler  Equation  of  Lateral  Motion 


V«5t  * 


v  4r  v«  *  244t  v;  • 24£v^  -  v;  *  *■ v 


P’(^/Tj) 


(34) 


0  <  %  <  1 
0  <  T| 


where. 


P  =  P.L-  ,  and 

K  El 

p*  =  p(^Ti) 


(35) 


and  its  boundary  conditions  are. 


v*  (0,Tj)  *  0  V^(1,T|)  =  0 

v*(0,T|)  =  0  v*u(l,Tp  =  0 


(36) 


Again,  discussion  of  the  initial  conditions  will  be 
delayed  until  later. 
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b.  Final  Euler  Equation  of  Axial  Motion 


L  +  (  ££  v-d^)  -  p  [-L  via,r\)  *  L  v;(l,Tl)]2  )  = 


1 

2pi; 


F(y\) 


(37) 


0  <  %  <  1 

0  <  T| 


and  initial  conditions, 


•MO)  =  L0 
i(0)  =  L0 


(38) 


F.  CASES 

There  are  two  general  cases  for  which  the  transient 
behavior  of  the  cantilever  beam  may  be  considered.  Recall 
that  the  beam  is  free  to  move  axially  when  an  axial  force  F(t) 
is  applied  resulting  in  an  axial  displacement,  or  when  an 
axial  displacement  L(t)  is  otherwise  imposed.  It  was  shown  in 
the  Euler  equation  of  axial  motion  (Eq.  37)  that  the  axial 
motion  described  by  L(t)  depends  on  F(t).  However,  if  this 
axial  motion  is  specified  simply  as  some  function  of  time 
alone  then  the  problem  is  greatly  simplified. 
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1.  Case  One,  L(t)  Prescribed 


If  L(t)  is  known  then  the  problem  is  reduced  to 
finding  a  solution  to  the  Euler  equation  of  lateral  motion 
(Eq.  34)  subject  to  its  boundary  and  initial  conditions.  The 
problem  is  a  linear,  initial-boundary  value  problem. 

An  even  further  simplification  occurs  if  the  axial 
motion  is  so  slow  that  the  time  derivatives  of  L(t)  are 
negligible.  Certain  linear  functions  of  L(t)  can  conveniently 
provide  such  a  condition  where  the  rate  (or  velocity)  is  made 
small,  and  the  second  time  derivative  (acceleration)  is  always 
zero.  In  this  case  Equation  (34)  reduces  to, 


v*  +  P^4  v*  -  0 
~eT  "  0 


,39) 


subject  to  its  boundary  and  initial  conditions  and  where 
p*(E,,rj)  =  0  for  a  beam  with  no  internal  loading.  Note  that  L4 
is  not  a  constant  here,  since  L(t)  is  a  prescribed  function  of 
time  which  needs  to  be  known. 

In  either  of  the  cases  (Eqs.  (35)  or  (39)),  a  closed 
form  solution  to  the  equation (s)  is  not  possible.  The  finite 
element  method  (FEM) ,  to  be  presented  in  Chapter  III,  was  used 
to  obtain  approximate  solutions  for  both  of  these  cases. 

2.  Case  Two,  F(t)  Prescribed 

In  the  case  where  F(t)  is  prescribed,  both  Equations 
(34)  and  (37) ,  subject  to  their  respective  boundary  and 
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initial  conditions  must  be  solved.  We  recall,  that  together 
these  two  equations  form  a  coupled,  initial-boundary  value 
problem.  Moreover,  they  are  both  now  nonlinear,  as  they 
contain  terms  with  both  dependent  variables,  L  and  D  and  their 
derivatives.  Therefore,  it  is  necessary  to  linearize  both 
equations . 

Any  number  of  different  strategies  are  possible  for 
the  linearization  of  these  equations.  The  strategy  used  here 
will  be  to  treat  each  dependent  variable  as  a  'primary'  or 
'secondary'  variable  in  accordance  with  the  following  scheme. 
The  assignment  of  'primary'  or  'secondary'  status  will  differ 
depending  upon  which  equation  is  to  be  linearized.  In  the 
linearization  of  the  equation  of  lateral  motion,  the  'primary' 
variables  are  lateral  deflection,  1)  and  its  derivatives;  and 
the  'secondary'  variables  are  axial  motion,  L  and  its 
derivatives.  As  'secondary'  variables  in  the  equation  for 
lateral  motion,  L  and  its  derivatives  are  evaluated  at  the 
previous  time.  On  the  other  hand,  for  the  equation  of  axial 
motion,  L  and  its  derivatives  are  considered  the  'primary' 
variables  and  V  and  its  derivatives  are  the  'secondary' 
variables.  In  this  case  V  and  its  derivatives  are  evaluated 
at  the  previous  time  step. 

The  linearization  of  Equation  (34)  is  quite  simple 
because  in  the  nonlinear  product  terms,  the  primary  variables 
(D  and  its  derivatives)  appear  linearly.  Therefore,  it  only 
becomes  necessary  to  evaluate  the  secondary  variables  in  these 
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product  terms  at  the  previous  time.  For  completeness  the 
linearized  equation  is  shown  below, 


(  .2  > 

*  P- 

$2 

i2 

L 

,  vti +  2C> 

i2 

~L 

Vl  ~  ~  Vl  +  L • 

V  J 

\  ) 

P* 

* 

(4/^1). 

0  <  ^  <  1 
0  <  T| 

(40) 

where  the  *  subscript  on  a  variable  (or  term)  denotes  that 
the  variable  (or  term)  is  evaluated  at  the  previous  time  and 
therefore  is  not  an  unknown  in  the  equation. 

The  linearization  of  Equation  (37)  is  not  so  simple 
because  the  primary  variables  (L  and  its  derivatives)  do  not 
appear  linearly  in  the  nonlinear  product  terms. 

If  Equation  (37)  is  expanded, 


L  + 


■zpi; 


V«(1,T1) 
-  - 


1_  [L2  v£2(1,T|)]  +  -L  [lL  v£  2  (1,T|)  •  V*  (1  ,T|)  ] 


■zz: 


(Tj) 


(41) 


each  of  the  bracketed  [  ]  terms  are  nonlinear.  These  terms 
can  be  linearized  in  a  number  of  different  ways.  Since  this 
is  primarily  an  "L"  equation,  the  "L"  operators  will  be 
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linearized  using  the  quasilinearization  technique  or  Bellman 


&  Kalaba  [Ref. 2].  The  nonlinear  terms  then  become, 


1. 


v«(l,T1> 

L2 


v«(  i,n) 


2 

where,  v^(l,T|)  = 

2.  L2  v^2(l,-n)  -  [0.(1, Ti)]2  {-l]  +  2L.l) 

3.  Li\\(  l,Tp  v'd,^)  =0(1,11)  V.(l,tl)  t.L 

4.  L2  V'2(1,T1)  »  (1,  TJ)  ( -L 2.  +  2  L.  L) 


~  V.(1,T1)  +  J.  e.d, ri) 

ll  2. 


(42) 


where  1.  is  the  length  of  a  finite  element  after  the  beam  is 
discretized,  0  represents  slope  \)^.  Recalling  that  t|=t,  in 
these  equations,  the  subscript  T|  denotes  partial 
differentiation  with  respect  to  time.  Again,  the  *  subscript 
on  a  variable  (or  term)  denotes  that  the  value  of  the  variable 
(or  term)  from  the  previous  integration  is  used. 
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Equation  (37)  can  be  rewritten  in  terms  of  the  L  operators 


and  the  following  groups  of  terms, 


Cx 

C2 

C* 


El  .  3 
~ZpL~o  L2. 


~ —2  V.  ( 1 ,  T|)  + 


El  .  2 

~ZpLo  IT 


V.  (1  ,TJ) 

• 


2  Ln 


e.2(i,ii)  i. 


-i-  e.(i,-n) 

+  J.  0.(1, TJ) 


c5  =  Jl  e. (i,rj)  v.<i,-n)  l. 

Lo 

C‘  =  -5X  v-2(1't') 
c7  =  -JL  V.2(1,T1>  L. 


CB 


1 

JpT0 


(43) 


Using  Equations  (43) ,  Equation  (37)  becomes, 

L  +  cx  +  C2L  +  C3  +  c4i  +  CSL  +  C6  +  C7L  =  CeF  (TJ)  (44) 

Equation  (40)  ,  its  boundary  and  initial  conditions, 
and  Equation  (44)  with  its  initial  conditions,  now  form  a 
coupled,  linear  initial-boundary  value  problem.  A  closed  form 
solution  of  these  equations  is  not  possible.  The  finite 
element  method  (FEM) ,  to  be  presented  in  Chapter  III,  could  be 
used  to  obtain  approximate  solutions  to  these  equations. 
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III.  SOLUTION  METHOD 


A.  FINITE  ELEMENT  METHOD  DEVELOPMENT 

Considering  the  L(t)  specified  case  first,  the  task  is  to 
solve  Equation  (34)  together  with  its  boundary  conditions 
given  by  Equation  (36) ,  and  its  initial  conditions  as 
determined  by  the  problem  being  investigated.  Definition  of 
the  initial  conditions  will  require  further  discussion  which 
will  be  conducted  later  in  this  development. 

Equation  (34)  is  a  linear  partial  differential  equation  in 
one  unknown,  v*(^,t)  when  L(t)  is  specified.  Recalling  that 
T|  =  t,  here  t  will  replace  TJ.  An  approximate  numerical 
solution  of  this  equation  together  with  its  initial  and 
boundary  conditions  can  be  obtained  by  a  Galerkin  finite 
element  formulation. 

1 .  Const  -auction  of  the  Beam  Element 

The  fourth  order  system  of  Equation  (34)  requires  C1 
continuity  (continuity  of  function  and  slope) .  In  order  to 
obtain  C1  continuity,  an  element  with  deflection  v',  and  slope 
0* ,  (0*  will  represent  in  the  development  that  follows;  as 
degrees  of  freedom  (DOF)  at  each  end  point  is  required.  This 
means  each  element  must  have  a  minimum  of  four  DOF,  which 
requires  a  cubic  polynomial.  These  interpolation  polynomials 
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are  the  set  of  cubic  spline  shape  functions  listed  below  and 
detailed  in  Appendix  A. 


,  3  ,  2  , 

q,  =  1  -  —2o?  *  a? 

q2  =  a  -  -La2  +  -La3 

■?. -  -  4«s 

J.  m 

q.  =  -_La2  +  _La3 

94  1.  ll 


(45) 


Figure  (A.l)  shows  that  shape  functions  qx  and  cq,,  are 
associated  with  the  displacement  DOF  (vj  ,\'2  where  subscripts 

1  and  2  represent  node  points  1  and  2)  at  the  element  end 
points;  and  the  even  numbered  shape  function  and  q4,  are 

associated  with  the  slope  DOF  (0i  ,  82)  at  the  same  locations. 

2  .  Global  FEM  Formulation 

In  terms  of  global  shape  functions  Q^,  the  FEM 
approximation  v*  to  V"  is  given  by, 

v*  *  v*  =  QT  8  =  £  Q,  5,  (46) 

i*l 


where  N  is  the  number  of  elements,  N  =  2 N  +  2  is  the  number  of 
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global  DOF,  and  6X  are  the  global  degrees  of  freedom.  The 
global  degrees  of  freedom,  8A  are  defined  as  follows, 

=  {8,8,  I  8,5.  i  5,8.  1  -  i  5Sl5ff>  (47) 

where  subscripts  1,2,3,  ..  .  (N- 1)  are  DOF  identifiers.  In  terms 
of  displacements  and  slopes,  we  have, 

5T  =  =  (vx  e,  :  v2  02  :  v3  e3  i  -  5  VH.X  0N.X  >*  (48) 

where  the  subscripts  1,2,3,  ...  (N+l)  refer  to  the  Global  Nodal 
Points  (GNP)  and,  <  >*  indicates  the  non  dimensional  forms  of 
and  0,  that  is  v*  and  0*. 

The  relationship  between  the  global  degrees  of  freedom 
5i  and,  v*  and  0*  ;  is  such  that  for  odd  i  (1 , 3, 5, . . .  N-l)  ; 

5i  =  9*  at  GNP  [  (i+1)  /2]  .  (i.e.,  §3  =  V*  (GNP  1),  53  =  (GNP 
2),  ...)  For  even  i  (2, 4,  6,  ...  N)  ;  =  0*  at  GNP  [i/2]. 

(i.e.,  52  =  0*  (GNP  1),  54  =  0*  (GNP  2),  ...) 

Each  ith  GNP  has  two  global  shape  functions,  and  hence 
two  DOF.  An  odd  numbered  shape  function  Qif  gives 

displacement  at  GNP{(i+l)/2]  =  vju.1)/2,  =  8i ,  and  an  even 
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numbered  shape  function  Qir  gives  slope  0*  at 
gnp  {i/2}  =  e;1/2,  =  di . 

3 .  Galerkin  FEM  Formulation 

In  accordance  with  the  Galerkin  Finite  Element  Method 
(FEM)  ,  we  form  the  approximate  solution  of  \), 

v*(4,t)  =  ?*(£, t)  =  QT(^)  5 (t)  (49) 

Using  the  above  approximation,  the  residual  function  for  the 
Euler  equation  of  lateral  motion  (Eq.  34)  is, 

R{k,t)  =  S£{9-}  -  P*  (50) 


or. 


(51) 


P 


2^-L  (V)4  -  2%t  (V*)%  -  ^L(V')^  +  LV* 
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After  the  final  substitution,  the  residual  is. 


P 


E2  *L  (qt  )  5 


R(k,t)  =  ( QT  )KK  6  + 

2^3  (<?t)4  5  -  2%L[q'  \  $ 
~  P' 


%L  [QT  \  6  +  L  QT  S 


(52) 


The  Galerkin  finite  element  equation  is  obtained  by  requiring 
that  the  residual  function  be  orthogonal  to  each  of  the  basis 
functions.  That  is. 


£  QR  d%  = 


(53) 


Substitution  of  Equation  (51)  into  Equation  (52)  gives, 

J> 0  (o' >«  •* 8  * 

£  %o(q’\ d%  5  -  2 p£  £  %ohe,\  at,  s  - 
Pi  £  %q  (qt  \  d%  s  *  p l  £  o  !ct  i  di,  8  =  £  o  p- 
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After  performing  two  integrations  by  parts  on  the  first  term 


(See  Appendix  D) ,  Equation  (53)  becomes, 

B. t.  .  /;  (o'  )«<*«*  P^  /;  5*  e  (o'  )K  d4  5  * 

2P^  £  %Q  (o'  \  5  -  2P£  £  %Q  (o'  \d\l  - 

P£  £  40  (o'  )t  d|  5  *  pi  £  0  (o'  )  d4  S  -  £  C  p-  dc, 


where  B.T.  is  the  vector  of  boundary  terms  resulting  from  the 
two  integrations  by  parts, 


B.T.  =  -  C>  (V*)^U 


(56) 


The  kronecker  delta  property  of  the  shape  functions  Qif 
results  in  the  reduction  of  the  boundary  term  vector  to, 


B.  T. 


-<**W°) 

(V’)K(0) 

0 

0 


0 

0 


■) 

(?*) 


(1) 

(D 


(57) 


31 


The  non  zero  terms  in  the  B.T.  vector  represent  the  non 
dimensionalized  loads  and  moments  applied  at  the  boundaries. 
If  there  is  no  applied  load  {or  moment)  at  the  boundary,  then 
B.T.  is  simply  the  zero  vector. 

A  discussion  of  the  second  term  of  Equation  (54) 
follows.  The  integration  by  parts  performed  on  the  first  term 
of  Equation  (54)  results  in  a  self  adjoint  (symmetric) 
operator,  a  condition  which  is  generally  desired  since  it 
reduces  storage  requirements  during  computer  processing. 
Integration  by  parts  on  the  second  term, 

JV  Q(QT)u  5  (58) 


results  in  a  B.T.  on  8*  evaluated  at  \  =  1.  Since  the  value 

of  8*  (1)  is  unknown,  an  integration  by  parts  was  not  performed 
on  the 
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Letting, 


a  -  /:  (o'  >«  < 

b  -  /;  r  0  (o'  >K  * 

c  =  £  ko{QT\  d%  (59) 

D  =  £  Q  (<?r  )4 


the  final  Galerkin  Equation  is, 


A  8  +  P—  B  8  +  C  8  -  <60> 

L  L 

2  pi  cS~Plc8  +  |3ld8  =  f 


The  details  of  the  construction  and  form  of  the  A,  B, 
C,  and  D  matrices  are  contained  in  Appendix  B. 

4 .  Integration  of  the  Galerkin  Equation 

a.  Reduction  of  the  Second  Order  System 

Equation  (59)  is  a  system  of  second  order  ordinary 
differential  equation  in  time.  For  numerical  integration 
purposes,  it  is  desirable  to  reduce  this  second  order  system 
to  a  first  order  system  in  time. 
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For  compactness,  let, 


JT=  -1 

1a  *  . 

f2£!  -  *] 

c 

L 

Lp  i 

- 

«=  ii 

L 


P  = 


(61) 


Then,  Equation  (59)  becomes, 


d5=m£+jc5+p 


(62) 


Letting  (0  =  & ,  it  follows  that, 


(63) 


and  Equation  (59)  now  becomes  the  first  order  system  of 
equations. 


h  =  © 

D(6  =  +  jc8  +  p 


(64) 
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In  explicit  matrix  form,  Equations  (62)  and  (63)  may  be 


written. 


1 

o 

M 

V 

1 

H 

o 

«1 

■> 

0 

1 

•  Q 

:  O 
_ i 

« 

CD 

II 

1 

X  i 

*  l 

1 _ 

< 

:  3: 

•  +  * 

p 

(65) 


Letting, 


I 

0 

r  *  =  ' 

M 

... 

•,  B  = 

1 

H  • 

o 

1 _ 

,  a.  =  • 

M 

... 

'  and,  £2  =  ' 

0 

0 

D 

w 

[K  i  M  ] 

W 

P 

(66) 


the  second  order  system  of  Equation  (59)  is  reduced  to  the 
following  first  order  system  in  time. 

=  (67) 


Vector  £2,  becomes  the  zero  vector  if  vector  P  is  a 
zero  vector.  Referring  to  Equations  (59)  and  (60)  we  see  that 
P  is  actually  defined  by  vector  F  which  is  defined  further  by 
the  boundary  term  vector  B.T.,  and  the  vector  describing  the 
contribution  of  an  internally  applied  load,  p*.  If  B.T.  is 
the  zero  vector  (no  applied  moment  or  load  at  the  boundaries) 


35 


and  there  is  no  internally  applied  load,  then  P  and  Q  are  zero 
vectors.  Evaluation  of  the  B.T.  vector  has  also  resulted  xu 
satisfying  the  natural  boundary  conditions  imposed  on  Equation 
(34) .  Equation  (66)  now  becomes, 


G  X  =  B  X 


(68) 


b.  Boundary  and  Initial  Conditions 

Prior  to  integrating  the  system  described  in 
Equation  (67),  the  boundary  and  initial  conditions  on  Equation 
(34)  must  be  imposed.  The  boundary  conditions  at  the  free  end 
(%=1) ,  were  imposed  through  the  boundary  term  vector  as 
previously  described.  The  strategy  used  to  impose  the 
essential  boundary  conditions  at  the  "fixed  end"  (^=0)  ,  is  one 
in  which  the  deflection  and  slope  at  the  "fixed  end"  are  set 
to  zero  when  the  X  vector  is  initialized,  and  the  G  and  H 
matrices  are  altered  such  that  the  deflection  and  slope  at  I; 
=  0,  remain  constant  with  time.  If  the  first  and  second  time 
derivatives  of  deflection  and  slope  at  £  =  0  are  constant  and 
equal  to  zero,  then  the  desired  conditions  of  zero  deflection 
and  slope  at  %  =  0  are  obtained  providing  that  the  initial 
conditions  on  deflection  (V(0,0)  =  0),  and  slope  (8(0,0)  =  0)  , 
are  satisfied. 

Initial  conditions  are  imposed  through  the 
initialization  of  the  X  vector  in  accordance  with  the  problem 
being  investigated.  Referring  back  to  the  global  FEM 
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formulation,  we  recall  that  each  nodal  point  has  two  DOF.  To 
satisfy  the  two  DOF,  deflection  (and  its  velocity)  ,  as  well  as 
slope  (and  its  velocity) ,  must  be  initially  defined  at  each 
.iodc .  The  iriitxal  cor.ditxCiis  also  satisfy  the  essential 
boundary  conditions  at  £  =  0 .  The  initial  conditions  are  more 
clearly  understood  if  the  X  vector  is  given  in  greater  detail, 


vi 

6! 


(69) 


where  the  *  subscript  indicates  the  terms  are  the  non 
dimensional  variables  (and  their  derivatives)  and  therefore, 
are  functions  of  t,,  not  x. 
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The  matrices  and  vectors  of  Equation  (67) ,  modified 


for  the  boundary  and  initial  conditions  follow, 


G  = 


i  1  0  0  -  0  0  0 
:  0  10-000 


(70) 


where  only  the  first  two  rows  of  D  are  altered  as  shown  and, 
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0  0  0  -  0  0  0 

0  0  0  -  0  0  0 


0  0  0  -  0  0  0 

0  0  0  -  0  0  0 


0  0  0  -  0  0  0 

0  0  0  -  0  0  0 


where  only  the  first  two  rows  of  I,  K,  and  M  have  been  altered 
as  shown  and. 
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and, 


X  = 


Vx 

VN- 

Vk 

Vi 

K 

V 

eK- 


(73) 


B .  FEM  VERIFICATION 

The  static  cantilever  beam  provides  a  problem  for  which  a 
known  solution  is  available  for  comparison  and  verification  of 
the  FEM  development  and  Fortran  code. 

For  the  static  cantilever  beam,  the  Euler  equation  of 
lateral  motion  (Eq.  1)  is  reduced  to, 

*JV„„(x)  =  p(x)  0  <x<L  (74) 
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with  the  boundary  conditions. 


V (0)  =  0 

vx(0)  =  0 


EJVXX(L)  =  M  (moment) 
=  P  (load) 


(75) 


Referring  to  the  Euler  equation  of  lateral  motion  (Fq. 
25) ,  and  considering  that  for  the  static  cantilever  beam, 


L  -  L  -  0  (76) 

8=0 


Equation  (73)  transformed  to  the  non  dimensional  coordinate  % 
becomes, 

=  0  <  k  <  i  <77> 


with  the  boundary  conditions, 

v (0)  =0 

v*<0)  =  0 


|fvu(D  =« 


(78) 


Referring  to  Equation  (34) ,  the  final  spatially  non 
dimensional  static  beam  equation,  where  the  lateral  deflection 
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has  also  been  non  dimensionalized,  becomes, 


Vtus<5>  -  -JjP<5) 


0  <  %  <  1 


(79) 


with  the  boundary  conditions, 


v*  (0)  =  0  ±tv^(l)  =  M 

Xi 

v*(0)  =  0  “v?K(l)  =  P 

L 


(80) 


The  Galerkin  FEM  formulation  for  Equation  (78)  and  its 
boundary  conditions  is  obtained  from  Equation  (59) .  Again, 
recalling  that, 


t  =  £  *  0  (81) 

8=0 


Equation  (59)  becomes, 


A  5  =  F 


(82) 


where, 


=  £  Op*  -  B.T. 


(83) 
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and,  if  there  is  no  excitation  internal  to  the  system  other 
than  at  the  boundaries  (p*  =  0)  then, 

F  =  0  -  B.T.  (84) 


The  boundary  conditions  given  in  Equations  (79) ,  must  be 
imposed  prior  to  solving  the  system  of  Equation  (81) .  To  do 
this  the  boundary  term  vector  B.T.,  resulting  from  the 

integration  by  parts  on  the  operator  and  shown  in 
Equation  (56)  , 


b.  r. 


-n«(0) 

v^(0) 

0 

0 

0 

0 

.-*«(D 


(56) 


is  evaluated  using  the  boundary  conditions  at  the  free  end  of 
the  beam  (£;  =  1)  . 
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Rearrangement  of  terms  in  Equations  (79)  gives, 


v«(D 


ML 

~EI 


v^d) 


(86) 


The  B.T.  vector  can  be  rewritten  as, 


B.  T. 


-r^{0) 

v^(0) 

0 

0 

0 

0 

PL2 
El 
_  ML 
El 

/ 


(87) 


Next  the  oundary  conditions  at  the  fixed  end  (£;  =  0)  are 
imposed.  Recalling  that  «  vj  =  v’(0)  and  52  «  0J  =  0'(O)  ,  the 

boundary  conditions  at  the  fixed  end  are  imposed  by  altering 
the  first  two  rows  of  both  the  A  matrix,  and  the  B.T.  vector 
to  force  and  82  to  zero. 
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The  final  system  can  be  written  as. 


10  0 
0  10 


N1 


0  0 
0  0 


(88) 


where  N  is  the  number  of  degrees  of  freedom. 
The  solution  to  this  system  is. 


6  = 


(89) 
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If 


PL2 

El 


(or  E) 
El 


is  set  at  unity,  then  to  obtain  the  8  for 


actual  values  of 


PL2 

~ei 


(or 


ML 

~EI 


the  8  vector  is  multiplied  by 


PL2 

'EY 


(or 


ML 

El 


) . 


The  exact  solutions  for  the  deflection  and  slope  at  the 
free  end  of  a  cantilever  beam  subject  to  a  concentrated  load 
(or  moment)  are  obtained  from  the  following  expressions, 


v  (L) 


_  PL3 
3 El 


e<L> 


_  PL * 


v  (L) 


ML2 

~ZEI 


e(L) 


ML 

~EI 


(90) 


where  P  (or  M)  is  the  load  (or  moment)  applied  at  the  free 
end. 

The  8  vector  is  the  vector  of  non  dimensional  deflections, 
and  slopes.  The  dimens  ial  vector  is  obtain  by  multiplying 
the  non  dimensional  displacements  (8*  for  odd  i)  by  L  in 
accordance  with  v  =  Lv*.  Since  slope  is  a  non  dimensional 
quantity  to  begin  with,  the  0*' s  (8A  for  even  i)  are  equal  to 
the  0's. 

The  Fortran  code  used  for  this  verification  and  comparison 
is  located  in  Appendix  C.  Results  of  the  verification,  also 
in  Appendix  C,  confirm  the  FEM  development. 
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C.  THE  F(t)  PRESCRIBED  SOLUTION  METHOD 

The  final  and  most  complex  case  posed  in  Chapter  II  was 
the  case  where  the  axial  force,  F (T)  is  prescribed.  In  this 
case  an  equation  for  axial  motion,  in  addition  to  the  Euler 
equation  of  lateral  motion,  was  required.  That  equation  was 
the  Euler  equation  of  axial  motion.  Together  these  two 
equations  form  a  nonlinear,  coupled,  initial-boundary  value 
problem.  After  the  linearization  of  these  equations,  the 
linear,  coupled,  initial-boundary  value  problem  consisted  of 
the  following  equations. 


P. 


f £2  A 


24 


v*  -  2%t.  -  %l,  v'  +  l,  v; 


■nn 


=  =  0 


(91) 


£  +  q  +  c2l  +  C3  +  C<£  +  CSL  +  C6  +  C7 L  =  CeFOD  (92) 

where  p*  =  0  for  the  no  internal  excitation  case,  and  Ci  are 
defined  in  Equations  (43)  . 
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By  defining  new  terms,  JA  which  require  updating  during 
the  time  integration  process, 

J,(t)  =  C2  +  C5  +  C7 

J2(t)  =  C4  (93) 

J 3<T)  -  Cx  +  C3  +  C« 

Equation  (45)  becomes, 

L  +  L  +  J2  i  =  CB  F(t)  -  J3  (94) 

Equation  (93)  is  a  second  order  differential  equation  in 
time.  Letting, 


S£  =  i 
a  =  L 


(95) 


and  using  substitution,  the  following  system  of  two  first 
order  differential  equations  is  obtain; 


a  =  L 

£  +  Ji  L  +  «J2  $£  =  C8  F(t)  -  J3 


(96) 


Since  L(t)  and  its  derivatives  do  not  depend  on  a  spatial 
variable,  Equations  (95)  do  not  require  a  FEM  formulation  and 
are  directly  added  to  the  system  of  equations  for  lateral 
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motion.  The  matrix  equations  for  lateral  motion  will  be 
similar  in  structure  to  the  system  of  Equation  (64) .  However, 
the  'secondary'  variables  (or  terms)  introduced  in  the 
linearization  of  the  equation,  are  assigned  their  values  from 
the  previous  integration.  Therefore,  the  sub  matrices  D,  K, 
and  M  which  contain  these  variables  (or  terms) ,  appear  with 
the  *  subscript.  £2  is  the  zero  vector  because  there  is  no 
internal  excitation,  and  the  natural  boundary  conditions  used 
to  evaluate  the  B.T.  vector  are  equal  to  zero  (moment =load=0)  . 
The  matrix  system  of  equations  which  is  obtained  after  adding 
Equations  (95)  to  the  system  in  Equation  (96)  is  shown  here. 


1 

o 

H 

i  ,,, 

[*1 

0  :  J 

[*1 

0 

. 

i 
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►  +  - 

... 
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o 

• 

1 _ 

CO 

K.  :  M, 

w 

0 

(97) 


49 


The  final  system  for  the  F(t)  prescribed  case  is, 


The  boundary  and  initial  conditions  for  the  equation  of 
lateral  motion,  and  the  initial  conditions  for  the  equation  of 


axial  mot ion ,  are  imposed  using  the  strategy  applied  in  the 
L(t)  specified  case,  and  described  in  Section  A  of  this 
Chapter . 
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IV.  CASE  STUDY  REPORT  AMD  CONCLUSIONS 


A.  GENERAL  DISCUSSION 

After  verification  of  the  finite  element  code  the 
investigation  of  the  transient  problem  began.  The  primary 
emphasis  of  the  studies  that  follow  is  on  obtaining  solutions 
to  problems,  and  not  on  investigation  of  numerical 
considerations.  However,  when  appropriate  the  researcher's 
thoughts  on  such  considerations  will  be  presented. 

The  case  studies  reported  are  investigations  of  the  L(t) 
prescribed  condition.  For  reference,  the  system  developed  in 
Chapter  III  for  the  L(t)  prescribed  case  is  repeated  here. 


f 

o 

H  : 

I _ 

< 

V 

►  — 

0  5  I 

4 

5 

u°  5  D 

(0 

K  l  M 

(100) 


or, 


Gi  =  HX 


(101) 


The  above  system  does  not  reflect  the  alterations  made  to 
impose  the  boundary  conditions  on  load  and  moment  at  the  free 
end  as  the  case  studies  addressed  only  the  case  of  homogeneous 
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boundary  conditions  (P=M=0) ,  with  no  internal  excitation 
(p*=0) .  Thus  the  P  vector  is  the  zero  vector,  and  does  not 
appear  in  the  3ystem  above. 

The  transient  problem  introduces  the  requirement  for  a 
numerical  integration  method.  To  perform  the  integration  on 
the  system  above,  the  IMSL,  Inc.  integration  subroutine, 
D I VP AG  was  chosen.  DIVPAG  is  a  double  precision  first-order, 
initial-value,  ordinary  differential  equation  solver. 

Two  classes  of  implicit  linear  multi  step  methods  are 
available.  The  first  is  the  Adams's  method  and  the  second  is 
the  backward  differentiation  formula  (up  to  fifth  order)  ,  also 
called  Gear's  stiff  method.  An  accepted  measure  of  stiffness 
is  the  ratio  of  the  maximum  and  minimum  eigenvalues  (X,„,  'Km±a) 
of  a  system.  A  problem  is  considered  stiff  for  very  large 
ratios.  The  vibrating  cantilever  beam  equation  of 
motion  results  in  a  stiff  system,  and  therefore,  Gear's  stiff 
method  is  used. 

1 .  Tima  Step  Convergence 

The  integration  routine  uses  an  internally  determined 
time  step  such  that  a  measure  of  global  error  does  not  exceed 
a  user  specified  tolerance.  This  feature  of  DIVPAG  provides 
error  control  to  the  user  of  the  integration  routine. 
However,  a  recognized  short  coming  of  this  integration  package 
as  applied  to  this  problem,  is  the  inability  to  update  the  G 
matrix  on  the  left  hand  side  of  Equation  (100)  at  each  of  the 
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subroutine  determined  time  steps.  That  is,  as  a  result  of 
the  call  structure  of  the  subroutine,  an  update  of  the  finite 
element  matrix  G  for  a  change  in  L(t)  and  its  derivatives,  is 
only  possible  outside  of  the  subroutine.  For  this  reason 
DIVPAG  is  placed  in  a  Do  Loop  and  the  G  matrix  is  updated  at 
each  entry  to  DIVPAG.  Although  the  H  matrix  could  be 
evaluated  inside  DIVPAG  via  a  FNC  subroutine  argument  of 
DIVPAG,  in  this  investigation  it  was  not.  It  was  updated  at 
the  same  time  the  G  matrix  was,  that  is,  at  each  entry  to 
DIVPAG.  The  accuracy  of  the  numerical  solution  depends  upon 
the  frequency  of  updating  of  the  G  and  H  matrices.  Entries  to 
DIVPAG  were  at  .025  second  intervals  for  all  case  studies  with 
the  exception  of  the  final  case  study  for  which  entries  were 
made  at  .01  second  intervals.  A  rapidly  changing  L(t) 
requires  more  updating  of  the  matrices  than  would  a  slowly 
changing  L(t).  In  effect,  a  solution  ultimately  should  be 
checked  for  "time  grid"  independence. 

2  .  Spatial  Grid  Convergence 

Convergence  of  the  solution  for  the  spatial  grid  is 
yet  another  consideration.  The  solution  is  a  function  of  the 
number  of  elements  (i.e.,  NDOF) .  For  linear  problems,  it  can 
be  shown  that  in  the  limit,  as  the  number  of  DOF  approaches 
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infinity,  the  approximate  solution  of  $  approaches  the  exact 
solution  V,  that  is, 


LIM  V  (102) 
NDOF  — >  °° 

However,  for  a  nonlinear  problem  there  is  no  guarantee,  but  a 
likelihood  that  the  approximate  solution  converges  to  the 
exact  solution  in  the  limit  as  the  number  DOF  approaches 
infinity.  A  preliminary  study  conducted  during  the  first  case 
study  showed  that  negligible  differences  existed  between  the 
eight  and  sixteen  element  solutions  for  that  particular 
problem.  This  was  the  basis  for  the  use  of  an  eight  element 
solution  for  all  subsequent  problems.  However,  it  is 
recognized  that  the  FEM  model  changes  with  length  (or  time) . 
Since  the  number  of  elements  (NEL)  is  constant  with  time, 
convergence  for  a  given  NEL  may  change  with  time  as  well. 
Furthermore,  the  effects  of  material  properties,  geometric 
dimensions,  and  functions  of  L(t)  (and  its  derivatives),  may 
also  influence  the  NEL  required  for  a  spatial,  grid 
independent  solution. 

3.  Computational  Effort 

Related  to  the  stiff  character  of  the  problem,  is  the 
very  large  amount  of  computational  effort  (CPU  time)  required 
to  obtain  solutions.  A  study  of  CPU  requirements  was  not 
conducted.  However,  integration  of  a  problem  over  a  real  time 
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five  second  period  took  as  long  as  a  week.  Typically,  DIVPAG 
performed  its  integration  over  2.N_6  second  steps.  Thus,  every 
.1  second  increment  in  time  required  approximately  500,000 
integration  steps.  Processing  was  conducted  using  the  Naval 
Postgraduate  School  IBM  3033  main  frame  time  share  system 
during  weekday  non  peak  hours  (1800-0900),  and  weekends. 
During  these  periods,  it  is  estimated  that  approximately  20 
percent  of  time  share  CPU  was  allocated  to  the  processing  of 
this  job.  A  restart  capability  was  coded  to  assist  in 
processing  during  non  peak  hours  only. 

4 .  The  Case  Study  Beam 

The  case  studies  that  follow  are  conducted  using 
material  properties  similar  to  those  of  plexiglass.  The 
modulus  of  elasticity  (E)  is  equal  to  100,000  psi.  The 
initial  length  of  the  beam,  L0  is  10.0  inches  for  all  case 
studies.  Two  moments  of  inertia  of  the  beam  cross-section 
(indication  of  beam  rigidity,  which  effects  the  stiffness  of 
the  problem)  are  obtained  from  the  cross-sectional  dimensions 
shown  in  Figure  (4.1).  In  recognition  of  the  stiff  nature  of 
the  problem,  and  in  interest  of  solving  a  realistic  problem  in 
the  minimum  amount  of  time  possible,  our  desire  was  to  select 
a  material  with  the  smallest  value  of  frequency,  which  is 
proportional  to. 


El 

jr* 


(103) 
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Thus,  for  a  beam  of  given  geometry  (I/L4)  ,  plexiglass  was 
selected  as  a  realistic  material  with  the  smallest  E/p  ratio. 


Figure  4.1  Beam  Geometry 

B.  FIRST  CASE  STUDY,  NEGLIGIBLE  DERIVATIVES  OF  L(t) 

The  strategy  used  in  the  case  studies  is  to  progress  from 
the  less  to  the  more  'difficult'  cases.  What  is  intended  by 
'difficult',  is  that  fewer  differential,  and  hence,  finite 
element  operators  are  involved  in  the  Fortran  coding  for  the 
less  difficult  case.  The  general  program  development  logic  is 
the  same  for  all  the  transient  case  studies,  however;  it  is 
generally  good  practice  to  limit  the  size  of  the  code  until 
the  logic  is  tested  and  functioning  as  expected.  By 
eliminating  the  derivatives  of  L(t),  only  the  A  part  of  the  K 
matrix  (See  Equation  (105))  and  the  D  matrix  of  Equation  (99) 
need  to  be  evaluated. 
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If  the  derivatives  of  the  specified  function  of  L(t)  are 


zero  or  negligible,  the  equation  of  lateral  motion  is, 


P£ l  v- 

EI  n 


=  0 


(104) 


subject  to  boundary  and  initial  conditions. 

For  this  case,  the  matrix  system  of  Equation  (99)  becomes, 


(105) 


where  the  K  sub  matrix, 


1 

k  =  ~-i 

1a  ♦  —B  * 

11!  -  £ 

C 

L 

L?  L 

L  J 

(106) 


is  reduced  to, 


K  = 


(107) 


This  case  was  examined  for  the  plexiglass  beam  with  the 
larger  cross  sectional  dimensions.  The  material  and  geometric 
properties  used  are  p (mass/unit  length)  =  6.988E-6  lbf’S2/in2 
(slugs/in)  and,  moment  of  inertia  I  =  81.38E-6  in4. 
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The  homogenous  boundary  conditions  are  imposed  as 
described  in  Chapter  III.  The  initial  conditions  (A.  vector) 
are  imposed  by  an  initial  parabolic  deformation  (Fig.  4.2)  of 
the  beam.  The  X  vector  is  initialized  for  all  DOF  according  to 
the  following  displacements  (and  slopes) , 


v*  (4* 0)  =  .1%2 

v;<s,o)  =  .24 


(108) 


and  velocities  of  displacements  and  slopes, 

V* (4, 0)  =  0  (109) 


Figure  4.2  Parabolic  and  Linear  Initial 
Conditions  Plot  Case  Study  One  and  Four 


The  prescribed  functions  of  L(t)  and  its  velocity  are. 


L(t)  =  L0  -  .  166fc 

£(t)  =  -.166 


(110) 


The  L (t ) ,  function  was  constructed  to  permit  the  beam  to  be 
drawn  half  way  (  lo  inches)  into  the  sleeve  (hole)  in  60 
seconds.  Figure  (4.3)  is  a  plot  of  L(t)  and  velocity,  £(t)  . 


time 


Figure  4.3  Length  and  Velocity  Function  Plot 

Case  Study  One 


This  problem  was  solved  for  a  four,  eight,  and  sixteen 
element  discretization.  This  was  the  only  investigation  of 
grid  independence  conducted.  The  results  of  this 
investigation  were  discussed  in  the  subsection  on  spatial  grid 
independence . 


C.  SECOND  CASE  STUDY,  PARABOLIC  FUNCTIONS  OF  L(t) 

Two  studies  are  conducted  where  L(t)  is  prescribed  by 
different  parabolic  functions.  If  L(t)  is  defined  as  a 
parabolic  function,  its  first  and  second  derivatives  are  non 
zero  and  significant  (See  Figures  4.5  and  4.6).  Thus,  the 
system  in  Equation  (99)  is  completely  defined.  In  addition  to 
the  dissimilar  functions  of  L(t),  the  two  cases  are  distinct 
in  their  cross-sectional  geometries. 

The  same  initial  conditions  were  imposed  for  the  two 
cases.  An  initial  deformation  of  the  beam  in  a  parabolic 
shape  was  imposed  again  as  in  the  first  case  study  and  again 
the  initial  velocities  are  zero.  The  A.  vector  of  Equation 
(99)  was  initialized  for  all  DOF  according  to  the  following 
displacements  and  slopes  (See  Figure  4.4), 


v* (£, 0)  =  %2 
vj(^,0)  =  2% 


(111) 
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(112) 


and  velocities, 

v*(4,0)  =  0 
v*(£,0)  =  0 

V,  V' 

2. 


1. 


0.0 

Figure  4  4  Parabolic  and  Linear  Initial  Conditions 

Plots  Case  Study  Two 

1.  Parabolic  L(t),  The  Less  Stiff  Beam 

"Case  One"  of  the  second  case  study  is  the  less 
"stiff"  problem.  Figure  (4.1)  shows  the  cross  sectional 
dimensions.  These  dimensions  result  in  material  and  geometric 
properties  such  that  p (mass/unit  length)  =  3.493E-6  lbf •  S2/in2, 
and  moment  of  inertia  I  =  10.17E-6  in4. 
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A  parabolic  function  of  L(t)  is  prescribed  such  that 
the  beam  is  drawn  to  half  its  original  length  in  2 . 5  seconds, 
reverses  direction,  and  returns  to  it  original  length  during 
the  next  2.5  seconds,  for  a  total  of  5  seconds.  Figure  (4.5) 
is  a  plot  of  L (t )  and  its  derivatives,  velocity  and 
acceleration.  The  functions  are, 

L(t)  =  L0  -  4t  +  .  8t2 

(113) 

L(t)  =  -4  +  l .  6t  {0  <  t  <  5.0)  sec. 

L(t)  =  1.6 


Figure  4.5  Length,  Velocity,  and  Acceleration  Plots 
Case  Study  Two  (Less  Stiff  Beam) 
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2.  Parabolic  L(t),  The  Stiff  Beam 

"Case  Two"  of  the  second  case  study  is  the  stiffer  of 
the  two  parabolic  L(t)  case  studies.  Figure  (4.1)  shows  the 
cross  sectional  dimensions  of  the  beam.  These  dimensions 
result  in  material  and  geometric  properties  such  that 
p  (mass/unit  length)  =  6.988E-6  lbf  •  S2/in2,  and  moment  of 
inertia  I  =  81.38E-6  in4. 

"Case  Two"  was  started  with  the  same  parabolic 
function  for  L(t)  as  "Case  One".  It  was  here  that  the 
significance  of  "stiffness"  and  computational  time  came  to 
focus.  Running  the  two  cases  simultaneously  as  separate  jobs 
on  different  system  accounts,  clearly  demonstrated  the 
difference  in  CPU  requirements  for  the  two  problems.  In  fact, 
there  was  such  a  disparity  in  computational  effort  that  it  was 
decided  to  change  the  course  of  the  stiffer  problem  ("Case 
Two")  such  that  it's  symmetric,  cycle  would  be  complete  in  2 . 1 
seconds  vice  the  5  seconds  of  "Case  One" .  The  functions  of 
L (t )  and  their  derivatives,  along  with  their  respective  time 
domains  are  given  here, 

L{t)  =  L0  -  4t  +  .  8t2 

L{t)  =  -4  +  1 . 6t  (0  <  t  £  1.0}  sec. 

Lit )  =  1.6 
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L(t)  -  33.2  -  50. 4t  +  24 

L  (t)  =  -50.4  +  48t  {1.0  £  t  £  1.1}  sec. 

L(t)  =  48.0 

L  (t)  =  5.128  +  .64  +  .  8t2 

£(t)  =  .64  +  1 . 6t  {1.1  <  t  <  2.1}  sec. 

L(t)  =  1.6 

These  functions  are  plotted  below  in  Figure  (4.6). 


Figure  4.6  Length,  Velocity,  and  Acceleration  Function  Plots 

Case  Study  Two  (Stiff  Beam) 
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3.  Discussion  of  the  Parabolic  Cases 


The  results  of  these  runs  provided  one  of  the  thought 
provoking  questions  of  the  research.  The  purpose  behind  the 
parabolic  function  of  L(t)  was  to  observe  the  beam's  activity 
in  the  case  where  it  returned  to  its  original  length  in  a 
symmetric,  cyclic  fashion.  The  interest  was  in  the  question, 
would  the  expected  symmetric  behavior  of  the  'pull-push' 
sequence  of  L(t)  be  predicted  by  the  code?  The  results 
clearly  showed  that  symmetry  did  not  occur  for  the  parabolic 
cases  (See  Figures  (4.10)  and  (4.11)).  In  fact,  the 
deflections  had  grown  considerably  during  the  'push'  stage  of 
the  L  (t )  cycle.  Where  did  the  energy  to  cause  such  large 
deflections  come  from?  One  possible  explanation  considered 
was  that  through  the  imposition  of  L(t),  energy  in  the  form  of 
work  had  been  added  to  the  system.  This  question  was 
addressed  in  the  final  case  study,  wherein  the  work  associated 

with  jFdL  was  tracked.  Another  possibility,  if  the  work 

cannot  account  for  the  increase  in  displacements,  is  that  the 
results  are  not  correct  due  to  a  break  down  in  the  numerical 
integration  during  the  latter  stage  of  the  'push'  stage  of  the 
cycle . 

Before  continuing  on  to  the  next  set  of  case  studies, 
a  discussion  of  what  is  a  most  thought  provoking  question 
resulting  from  the  research  thus  far  follows.  What  would 
happen  if  the  beam  where  drawn  totally  through  the 
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frictionless  hole  and  then  pushed  back  out  to  its  original 
length?  At  the  end  of  the  'puli'  stage  of  the  cycle,  the 
entire  beam  resides  motionless  in  a  straight  sleeve  (the 
frictionless  hole)  and  therefore  there  is  neither 
deformational  (strain)  energy  or  vibrating  (kinetic)  energy. 
Again,  energy  transfer  out  of  the  system  as  work  could  account 
for  this  phenomena.  In  any  case,  it  may  not  be  possible  to 
show  this  with  this  numerical  model  for  the  following  reasons. 
For  one  thing,  as  the  length  of  the  beam  shortens,  the  shear 
and  rotary  inertia  terms,  which  were  not  included  in  this 
model,  become  ever  increasingly  significant  and  in  fact  may 
dominate  the  physical  behavior.  Secondly,  even  if  the 
physical  model  could  be  modified  to  include  these  effects,  the 
frequencies  tend  toward  infinity  as  L(t)  approaches  zero,  and 
numerical  integration  would  not  be  possible. 

D.  HARMONIC  L(t)  PRESCRIBED 

Two  studies  were  conducted  simultaneously  on  two  beams 
with  the  material  and  geometric  properties  identical  to  the 
two  beams  used  in  the  parabolic  L(t)  study.  In  this  study, 
L  ( t )  was  prescribed  as  the  trigonometric  functions  of  L(t)  and 
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its  derivatives  given  here, 


Kt 

'1T3 


L(t)  =  9  +  COS  | 
i(t,  = 


L(t)  =  - 


7t 


1.5 


COS 


(^) 


{0  <  t  <  3.0}  sec.  (116) 


Accordingly,  L(t)  for  these  cases  varied  between  eight  and  ten 
inches  (Fig.  4.7).  The  symmetric,  cyclic  concept  was  used 
again  as  it  was  in  the  parabolic  L(t)  prescribed  cases. 


Figure  4.7  Length,  Velocity,  and  Acceleration  Function  Plots 

Case  Study  Three  (Both  Beams) 
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A  sinusoidal  function  was  also  prescribed  for  the  initial 
deformation  of  the  10  inch  plexiglass  beam.  The  X  vector  was 
initialized  for  all  DOF  in  accordance  with  the  following 
displacements  and  slopes. 


V(|,0) 

v;u,o) 


0.1 


1  ■  cos[£u 


7t 

20 . 0 


SIN\ 


<117) 


and  velocities, 


v* 0)  =  0  (118) 

vj<4»  0)  =  0 


The  initial  conditions  on  deflection  and  slope  are  shown 
graphically  in  Figure  (4.8) . 
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The  previous  results  of  increasing  displacements  (above 
the  initial  displacements)  observed  for  the  parabolic  L(t) 
studies  were  not  obtained  in  these  harmonic  L(t)  case  studies 
(See  Figures  (4.12)  and  (4.13)).  If  the  work  explanation  is 
the  correct  one  in  the  previous  section,  then  it  might  be  that 
work  is  associated  with  parabolic  L(t)  axial  motions,  and  not 
with  harmonic  L(t)  axial  motions.  In  order  to  investigate 
this  question  further,  a  investigation  was  undertaken  to  track 
work  for  the  parabolic  L(t)  case.  This  is  discussed  in  the 
next  section.  The  computational  effort  observations  of  the 
previous  cases  where  noted  again  as  well.  That  is,  the  CPU 
requirement  for  the  "stiffer"  problem  was  greater  than  for  the 
"less  stiff"  problem,  as  it  had  been  for  the  parabolic  L(t) 
cases . 

E.  CASE  STUDY  FOUR,  TRACKING  WORK  FOR  A  PARABOLIC  L(t) 

A  final  study  was  conducted  using  the  "less  stiff"  beam  in 
which  a  parabolic  L(t)  was  prescribed.  The  function  and  its 
derivatives  follow  and  are  plotted  in  Figure  (4.9). 

L(t)  =  10.0  -  2 . 66fit  +  .  886t2 

L(t)  =  2.666  +  l.lllt  (0  <  t  <  3.0}  sec. 

L  ( t)  =  1.771 
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The  A.  vector  was  initialized  for  all  DOF  according  to  the 
initial  deformation  of  the  beam  defined  by  the  following 
displacements  and  slopes, 


v*(4,0)  =  .l^2 
vj<$,0)  =  .2$ 


(120) 


and  velocities, 


v' (^, 0)  -  0  (121) 

vj<$,0)  =  0 


These  are  the  same  initial  conditions  as  used  in  the  first 
Case  Study  and  are  plotted  in  Figure  (4.2). 


Figure 


4 . 9 


time 


time 


time 


Length,  velocity,  and  Acceleration 
Case  Study  Four 


Plots 
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The  purpose  of  this  final  study  was  to  determine  whether 
the  increased  displacements  predicted  by  the  code  for  the 
parabolic  L(t)  cases-  during  the  latter  stages  of  the  'push' 
stage  of  the  cycle,  could  be  accounted  for  by  work  input  to 
the  system.  In  addition  to  tracking  work,  the  moment  and 
shear  were  also  tracked.  The  shear  diagrams,  shown  in  Figures 
(4.15)  to  (4.17),  and  moment  diagrams,  shown  in  Figures  (4.18) 
to  (4.20),  appear  to  be  reasonable.  The  small  values  of  these 
parameters  is  due  to  the  values  of  Young's  modulus,  E,  and 
moment  of  inertia,  I,  used  in  this  study.  The  product  of  El 
for  the  cases  studies  here  are  1.017  lb.  in2,  and  8.138  lb. 
in2 . 

The  diagrams  for  axial  force  F  and  work  W,  shown  in  Figure 
(4.21)  do  not  appear  to  be  reasonable  and  therefore  are 
suspect.  Assuming  a  one  to  one  relation  between  L(t)  and  F(t) 
exists,  it  is  difficult  to  imagine  that  such  a  force  would 
produce  the  smooth  parabolic  L(t)  and  vice  versa.  A  tentative 
conclusion  therefore  is  that  either  F  (t)  was  not  coded 
correctly  or  that  there  is  in  fact  an  instability  in  the 
numerical  integration  during  the  latter  stage  of  the  'push' 
cycle  of  the  problem.  An  effort  is  presently  underway  to 
determine  if  the  coding  for  the  calculation  of  F(t)  is 
correct.  It  should  be  remarked  however  that  prior  to  the 
erratic  behavior  of  F(t),  which  occurs  late  in  the  'pull-push' 
cycle,  the  values  of  F(t)  seemed  reasonable. 
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F.  FINAL  COMMENTS  AND  RECOMMENDATIONS 


The  results  of  this  initial  investigation  on  the  behavior 
of  a  vibrating  beam  subject  to  a  prescribed  axial  motion  leads 
to  the  following  conclusions.  First  and  foremost  is  that  the 
implementation  of  the  FEM  numerical  scheme  was  accomplished 
with  success,  although  it  is  not  certain  that  some  numerical 
difficulties  are  not  encountered  at  the  later  stages  of  the 
analysis.  Further  work  must  be  undertaken  to  resolve  whether 
the  increase  in  vibration  amplitude  which  is  predicted  by  the 
code  is  an  actual  result  of  work  input  to  the  system  or 
whether  it  is  associated  with  a  numerical  instability.  Prior 
to  the  investigation  of  the  'real'  gun  barrel  problem,  One 
might  also  investigate  whether  the  omission  of  axial  strain 
energy  form  the  model,  which  is  common  whenever  bending  and 
bar  activity  coexists,  could  account  for  this  behavior.  If 
so,  additional  terms  for  axial  strain  energy  could  be  included 
in  the  formulation. 

It  is  interesting  to  note  that  the  equation  of  axial 
motion  relates  the  axial  force  F(t)  not  only  to  the  axial 
acceleration  L,  in  accordance  with  Newton's  law  of  motion  for 
rigid  bodies,  but  also  includes  additional  terms  associated 
with  the  deformational  strain  energy  of  bending,  and  the 
kinetic  energy  of  beam  vibration,  at  the  free  end  of  the  beam. 
The  former  term  adds  to  the  acceleration  term  while  the  latter 
term  decreases  it. 


73 


The  practical  problem  associated  with  the  axial  motion  of 
a  gun  barrel  due  to  the  recoil  action  of  firing,  which 
provided  the  impetus  for  this  study,  was  formulated  but  not 
solved  here.  An  experimental  investigation  should  be 
undertaken  to  ascertain  the  accuracy  of  the  numerical  model. 
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PARABOLIC  AXIAL  MOTION 


Figure  4.10  Parabolic  Axial  Motion  Transient  Response 
Case  Study  Two  (Stiff  Beam) 
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PARABOLIC  AXIAL  MOTION 


Figure  4.11  Parabolic  Axial  Motion  Transient  Response 
Case  Study  Two  (Less  Stiff  Beam) 
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HARMONIC  AXIAL  MOTION 


PARABOLIC  AXIAL  MOTION 


TIME  (SEC) 


Figure  4.14  Parabolic  Axial  Motion  Transient  Response 

Case  Study  Four 


79 


APPENDIX  A 


THE  CUBIC  SPLINE  SHAPE  FUNCTIONS 

The  beam  element  is  constructed  using  four  shape  functions 
(qlf  CJ2 f  ^3/  and  ^4)  which  satisfy  the  following  conditions. 

q1(NPl)=l  ,  q[  (NP1 )  ~  0 

g1(NP2)=0  ,  q'1(NP2)=0 

I  2. 

q2  (npi)  =  0  ,  g2  (wpi )  - 1 

q2  (NP2)  -  0  ,  q2(NP2)=0 

q3  {NP1)  =  0  f  q3  {NP1 )  =  0 

J:  g3  ( NP2 )  =  1  ,  g£(MP2)  =  0 

»  Z 

?4  (NPI)  =  0  ,  g«  (NPl )  =  0 

g4(WP2)=0  ,  g*(WP2)=l 

where  the  (')  superscript  represents  a  differentiation  with 
ieSj^ect  to  the  spatial  variable. 
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By  satisfying  the  four  conditions  on  each  q±,  the  four 
element  shape  functions  can  be  constructed  from  a  cubic 


equation  of  the  form, 

gi  =  ai  +  bi  a  +  ci  a2  +  d±  a3 

qx  is  constructed  here  as  an  example. 

q1  =  ax  +  bx  a  +  Cj  a2  +  dx  a3 

1.  q1  (NP1)  =1  =>  1  =  a1 

2.  q[  (NP1)  =0  =»  0  =  b, 

ii 

gx  =  1  cx  a2  +  dx  a3 

U 

3.  q1  (NP2)  =0  =*  0  =  1  +  cx  l2m  +  dx  l\ 

4.  qi  (NP2)  =  0  =>  0  =  2cx  1,  +  3di  l\  =$  cx  =  - 

.  2  3 

substitute  (4.)  — >  (3.)  =>  dx  =  — _  and,  cx  =  - 

•  -Z# 

thus, 

gx  =  1  -  ^  a2  +  a3 
^  ll 
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n  (cm 


Using  the  conditions  which  define  shape  functions  42,  q3, 
and,  q4  (listed  on  the  previous  page)  the  other  three  shape 
functions  are  obtained. 
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APPENDIX  B 


CONSTRUCTION  OF  THE  GLOBAL  MATRICES 


The  global  matrices  A,  B,  C,  and,  D  are  constructed  from 
the  element  matrices  a*,  b®,  c®,  and,  d®  according  to 

relationships ; 


A  =  £  Q"{QT  )"d%  =  U  a® 

B  =  £  Z,2  Q  {qt  )"dE,  =  Ub® 
C  =  f1  h,Q  ( QT  )'<*%  =  U  c® 

Jo 

d  =  f  Q  (QT  )' d%  =  U  d® 

Jo 


(129) 


where, 

a®  =  £  q"{qT  )"  da 
b .  =Vjz  q  ( qT  )"  da 

(130) 

c-  =  ljx  q  ( qT  )'  da 
d*  =  q  ( qT  )  da 


given  that  E,  is  a  constant  which  approximates  E,  transformed 
to  the  local  coordinate  a. 
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The  transformation  of  H,  is  as  follows.  Referring  to 
Figure  (B . 1 ) , 

-,1 

^  '  h.  where,  =  £  1 

i-  1  ' 

and,  H,2  =  a2  •  2 a  +  r|* 

Because  the  transformation  of  £,  to  a  function  of  tt  results  in 
integrals  which  are  difficult  to  evaluate,  it  is  desireable  to 
use  an  alternative  strategy.  If  we  let, 

%  -  l  =  M,  *  4= 

a  numeric  value  can  be  assigned  to  this  quantity.  Thus,  the 
difficult  integration  is  eliminated.  Any  accuracy  lost  in  the 
approximation  will  be  recovered  by  additional  interations  to 
obtain  convergence  of  the  FEM  solution. 


t- 

•- 


♦ 


1  '<4 


♦ - \) - * - • 


F igure  B . 1 
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The  final  a*,  b",  c*,  and,  d*  matrices  follow, 


6 

_  11 

6 

1 

* 

-H 

IT) 

10 

51. 

To 

1 

_21- 

1 

1. 

'T 3 

15 

IT 

6 

1 

6 

11 

51. 

To 

51. 

10 

1 

1. 

1 

• 

CM 

TT 

FT 

TF 

TF 

92 


1 

l. 

1 

±. 

T 

IT 

T 

IT 

1. 

0 

TT 

Tft 

"ST 

l 

l 

■i. 

~2 

TT 

~1 

TT 

l 

ll 

l 

• 

TO 

0 

131. 

m2 

91. 

1312 

35 

210 

~w 

420 

lll2 

i2 

13  ll 

_  i2 

TlTT 

TT5 

TST 

140 

13 12. 

131. 

m2 

“7TT 

“33" 

~2TTT 

131^ 

_  i2 

_iii: 

i2 

TTT 

T4T 

"Tnr 

T55 
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APPENDIX  C 


STATIC  CANTILEVER  BEAM  FORTRAN  CODE 
AND  SAMPLE  OUTPUT 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX X 
x  MARK  R.  DEVRIES  IT  USCG 

x  NAVAL  POSTGRADUATE  SCHOOL 

x  SEPTEMBER  1990 

x  THESIS 

x  MASTER  OF  SCIENCE  IN  MECHANICL  ENGINEERING 
x 

x  TITLE: 

X  VIBRATION  OF  A  CANTILEVER  BEAM 

x  THAT  SLIDES  AXIALLY  IN  A  FRICTIONLESS  HOLE 
x 

x  THE  FOLLOWING  FORTRAN  CODE  IN  A  VERIFICATION  OF  THE  FINITE 
x  ELEMENT  FORMULATION  FOR  THE  TRANSIENT  PROBLEM  TO  BE 
x  PURSUED  IN  THE  NEXT  PROGRAMING  STEP, 
x  THE  PROGRAM  VERIFIES  THE  INITE  ELEMENT  METHOD 
X  CODE  LOGIC  ON  THREE  POSSIBLE  STATIC  BEAM  PROBLEM. 

X  (1)  FIXED  END  WITH  TWO  ROLLER  SUPPORTS,  ONE  AT  THE  CENTER 
X  AND  A  SECOND  AT  THE  OPPOSITE  END.  THIS  BEAM  IS  LOADED  BY 
x  A  CONCENTRATED  MOMENT  AT  THE  ROLLER  SUPPORTED  END. 

X  (2)  AND  (3)  ARE  CANTILEVER  BEAM  PROBLEMS,  ONE  LOADED  BY  A 
X  CONCENTRATED  LOAD  AT  THE  FREE  END  AND  THE  OTHER  LOADED 
X  BY  A  CONCENTRATED  MOMENT. 

x  THE  PROGRAM  IS  NOT  FLEXIBLE  IN  THAT  IT  REQUIRES  EDITTING 
X  AS  NOTED  IN  COMMENT  LINES  IN  THE  FOLLOWING  SUBROUTINES 
X  DEPENDING  ON  WHICH  OF  THE  (3)  CASES  IS  BEING  RUN.  THE  SOLE 
X  PURPOSE  OF  THIS  PROGRAM  IS  TO  VERIFY  THE  FEM  FORMULATION, 
x  IT  IS  NOT  INTENDED  TO  IMPRESS  SOLFWARE  ENGINEERS, 
x 

X  (1)  SUBROUTINE  BC 
X  C  2 )  SUBROUTINE  OUTPUT 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

x  VARIABLE  IDENTIFICATION  X 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

x  NEL  -  NUMBER  OF  ELEMENTS 

x  N5.NP  -  NUMBER  OF  SYSTEM  NODAL  POINTS 

x  NDOF  -  NUMBER  OF  DEGREES  OF  FREEDOM 

x  E  -  MATERIAL  MODULUS  OF  ELASTICITY 

x  GI  -  SECOND  MOMENT  OF  THE  BEAM  CROSS-SECTION  AREA 

X  BLGTH  -  BEAM  STATIC  LENGTH 

x  ELE  -  ELEMENT  LENGTH 

X  BCM  -  EXTERNALLY  APPLIED  MOMENT  AT  FREE  OR  SIMPLY  SUPPORTED 
x  END 

x  BCFORC-  EXTERNALLY  AFPLIED  FORCE  AT  FREE  OR  SIMPLY 
x  SUPPORTED  END 

x  NDETRM  -  VARIABLE  USED  IN  LOGIC  STATEMENT  FOR  TYPE  OF  B.C. 
x  SLOPE  -  SLOPE  AT  FREE  END  OF  CANTILEAVER  BEAM 
X  DEFLEC  -  DEFLECTION  AT  FREE  END  OF  CANTILEAVER  BEAM 
x  FACTOR  -  SCALOR  NON-DIMENSIONAL  GROUP 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


PARAMETER  (N=70) 


ooono  oooo  noon 


DIMENSION  A(N,N) ,F(N) 

DIMENSION  WKAREA(6000) 

OPENtlO,  FILE  =  '/MATRIX  OUTPUT') 

CALL  ZERO  r A , F) 

CALL  DATA  ( NEL , NSNP , NDOF , E , GI , BLGTH , EL E, BCM, BCFORC ) 

CALL  MATX  ( A , F , ELE, NEL , NDOF) 

CALL  BC  (A,F,NDOF,NDETRM) 

'■IRITE  (10,50  ((A(J,K),K=1, NDOF), J=l, NDOF) 

CALL  LE0T2FCA, 1 , NDOF, N, F , 0 , WKAREA, IER ) 

CALL  OUTPUT  ( NSNP, ELE, F, BCM, BLGTH, E, GI , NDETRM, NDOF, BCFORC) 

STOP 

END 


ZERO  ALL  MATRICES 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
SUBROUTINE  ZERO  ( A , F ) 

PARAMETER  (1070) 

DIMENSION  A ( N , N ) ,  F(N) 

DO  20  1  =  1,  N 
F( I )  =0 . 0 
DO  10  J  =  1 ,  N 
A( I , J  )  =  0 . 0 

1C  CONTINUE 

20  CONTINUE 

RETURN 
END 

INPUT  DATA 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
SUBROUTINE  DATA  ( NEL , NSNP, NDOF, E, GI , BLGTH, EL E , BCM, BCFORC) 

PRINT  x,  'ENTER  THE  NO.  OF  ELEMENTS  TO  BE  USED  IN  THE  APPROX 
READ  x , NEL 
WRITE! 6 , 20 )  NEL 

20  FORMAT  (/2X,'N0.  OF  ELEMENTS  IS', 15) 

NSNP=NEL+1 
WRIT  £( 6 , 25 1  NSNP 

25  FORMAT ( /2X , ' NO .  OF  SYSTEM  NODAL  POINTS  IS', 15) 

ND0F=2XNSNP 

WRITE! 6,76)  NDOF 

26  FORMAT ( / 2X , ' NO .  OF  D.O.F.  IS', 15) 

PRINT  x , 'THE  MODULUS  OF  ELASTICITY  IS?' 

READ  X , E 

WR I T E( 6 , 27  )  E 

27  FORMAT  (/2X, 'MODULUS  OF  ELASTICITY  IS',  F10.1) 

PRINT  X.'THE  SECOND  MOMENT  OF  THE  BEAM  CROSS-SECTION  AREA  IS 

READ  x , GI 

WR I TE( 6 , 28  )  GI 

23  FORMAT ( / 2X , 'THE  SECOND  MOMENT  IS',F10.1) 

PRINT  x , 'THE  INITIAL  LENGTH  OF  THE  STATIC  BEAM  IS?' 

READ  x,  BLGTH 
WRITE! 6 , 29  )  BLGTH 

29  F0RMAT(/2X, 'THE  BEAM  LENGTH  IS',F8.3) 

EL  E  =  1 . O/FLOAT ( NEL ) 

PRINT  x, 'ENTER  THE  VALUE  OF  THE  APPLIED  MOMENT' 

READ  x , BCM 
WRITE  (6,30)  BCM 

30  FORMAT (/2X, • MOMENT =' , F8 . 1 ) 

PRINT  x, 'ENTER  THE  VALUE  OF  THE  APPLIED  FORCE* 

READ  X , BCFCRC 
WRITE  (6,40)  BCFORC 
40  FORMAT (/2X, ’ FORCE= ' , F8 . 1 ) 

RETURN 

END 

FILL  LARGE  A 

XXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

SUBROUTINE  MATX  ( A  F, EL E, NEL , NDOF ) 

PARAMETER  (N=70) 
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DIMENSION  A(N,N),F(N) 

C 

C  CALCULATE  LITTLE  A  MATRIX 
C 

A11=12.0/CELE**3.0) 

A12=6.0/(ELE*#2.0) 

A13= ( -1 .0 )  XA11 
A14=A12 
A21=A12 
A22=4. O/ELE 
A23=(-l  .  0)xA12 
A24  =  2 . O/ELE 
A3I =  A1 3 
A32=A23 
A33=A11 
A34=A23 
A41  =A14 
A42= A24 
AA3-A3A 
A44=A22 
C 

C  FILL  LARGE  A  MATRIX 

C 

L  =  2*NEL-1 
C 

DO  10  1=1, L, 2 
C 

AC  I, I )=A( I , I )  +  All 
AC  1 ,  1  +  1  )=A( 1 , 1  +  1  )+A12 
A  C I , I  +  2)  =  A13 
AC  I , I+3)=A14 
C 

ACI+1,I)=ACI+1,I)+A21 
ACI+1,I+1)=ACI+1,I+1)+A22 
ACI+l, I+2)=A23 
AC  1  +  1, I+3)=A24 
C 

AC  1  +  2 , 1 ) =A31 
AC  1  +  2, 1  +  1  )=A32 
ACI+2,I+2)=A33 
A ( I +2 , I +3 ) =A34 
C 

AC  1  +  3, 1 )  =  AA1 
ACI  +  3, 1+1  )=AA2 
ACI+3, I+2)=AA3 
AC  1  +  3, I  +  3)=AAA 
C 

10  CONTINUE 
C 

RETURN 

END 

C 

C 

C  *#*##¥*X**X¥¥*X**X**X#*#»#XXX##**XX#«#«X#*X#*#*X*#****#**XX# 

C  THE  FOLLOWING  SUBROUTINE  ALTERS  THE  GLOBAL  MATRICES  TO 

C  IMPOSE  THE  BOUNDARY  CONDITIONS 
C 

SUBROUTINE  BC  C A , F , NDOF , NDETRM) 

P ARAMET ER C  N  =  7  0 ) 

DIMENSION  A ( N , N ) , F C  N ) 

PRINT  x, 'ENTER  1  FOR  THE  OVER  DETERMINANT  CASE  OR' 

PRINT  x; ' 2  FOR  THE  FREE  END  CASE.' 

READ  x,  NDETRM 
IF  C  NDETRM  .NE.  1)  GOTO  20 
C 

C  AMEND  A  TO  ACCOUNT  FOR  BOUNDARY  CONDITIONS 

C 

C  CHANGE  FIRST  AND  SECOND  ROWS  TO  ACCOUNT  FOR  THE  ESSENTIAL 
C  BOUNDARY  CONDITIONS  AT  THE  FIXED  END. 

C 

C 

C  THE  J-TH  EQUATION  IS  THE  EQUATION  DESCRIBING  DEFLECTION  AT 
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nnnnnnnnnn  on  oooooooo  noooooo  oonoo  o  o  o  o  ooooooo 


THE  LOCATION  OF  THE  CENTER  SUPPORT .  IT  IS  REPLACED  BY  THE 
ESSENTIAL  BOUNDARY  CONDITION  ON  DEFLECTION. 

THE  (NDOF-l)TH  EQUATION  IS  THE  EQUATION  DESCRIBING  THE 
DEFLECTION  AT  THE  ROLLER  SUPPORTED  END.  THIS  EQUATION 
IS  REFUGED  BY  THE  ESSENTIAL  B.C.  ON  DEFLECTION. 

J=ND0F/2 

DO  10  I = 1 , NDOF , 1 

A( 1 , 1 )  =  0 . 0 
A  (  2 ,  I )  =  0 . 0 
A  (  J  ,  I )  =  0 . 0 
A( NDOF-1 ,  I )  =  0 . 0 

10  CONTINUE 

A(l,l)=1.0 
AC2,2)=1.0 
ACJ, J)=l  .0 

A(N D0F-1, ND0F-1)=1  .0 
GO  TO  40 

20  DO  30  1=1, N DCF, 1 
A( 1 , I ) =0 . 0 
A  C  2  ,  I  )  =  0  .  0 
30  CONTINUE 
A  C 1 , 1 )  =  1  .0 
A ( 2 , 2 )  =1 . 0 

THIS  LINE  IS  ACTIVATED  FOR  THE  CANTILEVER  BEAM  LOADED 
BY  A  CONCENTRATED  LOAD  CASE 

F ( NDOF-1 ) = 1 . 0 

xxxxPOINT  MOMENT**#* 

THIS  LINE  MUST  BE  ACTIVATED  FOR  BOTH  THE  OVER  DETERMINATE  CASE 
AND  THE  CANTILEVER  BEAM  WITH  A  CONCENTRATED  MOMENT  CASE 

F( NDOF) =1 .0 

40  RETURN 
END 


FORMULATE  OUTPUT 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
THIS  SUBROUTINE  CALCULATES  THE  EXACT  SOLUTION  FOR  THE 
CANTILEVER  BEAM  CASES  AS  WELL  AS  PRINTS  THE  OUTPUT  OF  ALL 
CASES. 

SUBROUTINE  OUTPUT  ( NSNP , ELE, F, BCM, BLGTH, E> GI , NDETRM, NDOF, BCFORC) 
PARAMETERC  N  =  7  0 ) 

DIMENSION  FCN) 

xxxPOINT  FORCE#*# 

FACT0R=BCF0RCX( BLGTHXX2 . 0)/( E#GI) 

DO  5  1=1, NDOF-1, 2 

F  C I )  =  F ( I )*BLGTH#FACTOR 
F(I+1)=F(I+1)*FACT0R 
5  CONTINUE 

xxxPOINT  MOMENT*** 

FACT OR =  BLGTH* BCM/ ( E#GI ) 

DO  5  1=1, NDOF-1, 2 

F( I )=F( I ) *BLGTH*FACTOR 
F( 1+1 ) =F( 1+1 )* FACTOR 
5  CONTINUE 

MRITE( x , 30  ) 

WRITEL  # , 40  ) 
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onnoooo  oooo 


IF  ( NDETRM  .NE.  1)  GOTO  15 
HRITEC  *, 50) 

J  =  1 

XL  0C=  0 . 0 

DO  10  I=1,NSNP,1 

WRIT  EC  * ,  20 )  XLOC,  FCJ),  FCJ+l) 

XL0C=XL0C+ELE 
J  =  J  +  2 
10  CONTINUE 
GO  TO  80 

#xx*POINT  LOAD**#* 

ACTIVATE  FOR  THE  CANTILEVER-CONCENTRATE  FORCE  CASE 

15  SLOPE=CBCFORC*(BLGTH**2.0) )/C2. 0*E*GI  ) 
DEFLEC=(BCFORC*(BLGTH#*3.0))/(3.0#£*GI) 

****POINT  MOMENT***** 

ACTIVATE  FOR  THE  CANTILEVER  BEAM  CONCENTRATED  MOMENT  CASE 

15  SLOPE-(BCM*BLGTH)/(E#GI) 

DEFLEC=CBCM*( BLGTH**2 . 0 ) )/ ( 2 . 0*E*GI ) 

HR I TEC #,60) 

HR I T  E ( * , 6  5 ) 

WRIT  EC  * , 7  0  )  DEFLEC, FCNDOF-l), SLOPE, F(NDOF) 

FORMAT (2X, F8. 3, 3X,E12. 4, 3X,E12. 4) 

FORMAT <2X, »  X-LOCAT ION • , 3X , 'DEFLECTION' , 3X, 'SLOPE* ) 

FORMAT ( IX, '  DEFLECTION  AT  B  SLOPE  AT  B') 

FORMAT (6X,  'EXACT MIX,  'FEMM1X,  'EXACT' ,  11X,  'FEM* ) 

FORMAT (1X,E12.4,3X,E12.A,3X,E12.4,3X,E12.A) 

RETURN 
END 


20 

50 

60 

65 

70 

SO 
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NO.  OF  ELEMENTS  IS  8 

NO.  OF  SYSTEM  NODAL  POINTS  IS  9 

NO.  OF  D.O.F.  IS  18 

MODULUS  OF  ELASTICITY  IS  30000000.0 

THE  SECOND  MOMENT  IS  100.0 

THE  BEAM  LENGTH  IS  100.000 

ENTER  THE  VALUE  OF  THE  APPLIED  MOMENT 

MOMENT  =  0.0 

ENTER  THE  VALUE  OF  THE  APPLIED  FORCE 
FORCE=  1000.0 

ENTER  1  FOR  THE  OVER  DETERMINANT  CASE  OR 
2  FOR  THE  FREE  END  CASE. 

? 

2 

DEFLECTION  AT  B  SLOPE  AT  B 

EXACT  FEM  EXACT  FEM 

0.1111E+00  0.1111E+00  0.1667E-02  0.1667E-02 
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APPENDIX  D 


INTEGRATION  BY  PARTS 


The  following  is  the  detail  of  the  integration  by  parts  on 
the  first  term  of  Equation  (54) , 


/:  o  u  « 8 


(137) 


The  first  integration  results  in, 

/:  o  (»'  u  « 6  - 

(138) 

0  (O'  W  8  -  £  0t  (O'  )w  *  » 


A  second  integration  performed  on  the  integral  in  Equation  (2) 
gives. 


-  r  o,  <o'  )m  <« « - 

-  0,  (O'  j;  0„  (O’  )«  4,  8 


(139) 
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Combining  Equations  (1)  ,  (2)  ,  and  (3)  gives  the  symmetric 

operator  and  boundary  terms  below, 


/;  0  (C'  )UK  *  >  - 

9  (°r  L  5  -  <eT  >«  6 1  *  «« <o'  >«  8 


(140) 


0 
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APPENDIX  E 


TRANSIENT  BEHAVIOR  OF  A  CANTILEVER  BEAM  FORTRAN  CODE 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxx 
X  MARK  R  DEVRIES  IT  USCG 

X  NAVAL  POSTGRADUATE  SCHOOL 

X  SEPTEMBER  1990 

X  THESIS 

X  MASTER  OF  SCIENCE  IN  MECHANICAL  ENGINEERING 
x 

X  TITLE; 

X  VIBRATION  OF  A  CANTILEVER  BEAM 

X  THAT  SLIDES  AXIALLY  IN  A  FRICTIONLESS  HOLE 
x 

x  THE  FOLLOWING  FORTRAN  CODE  UTILIZES  THE  FINITE  ELEMENT 
X  METHOD  AND  AN  IMSL  PACKAGE  INTEGRATION  SUBROUTINE  DIVPAG 
X  TO  SOLVE  THE  ABOVE  PROBLEM.  THE  PROGRAM  IS  WRITTEN  WITH 
x  NUMEROUS  COMMENT  LINES  WHICH  EXPLAIN  THE  CODING, 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  VARIABLE  IDENTIFICATION  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  NEL  -  NUMBER  OF  ELEMENTS 
X  NSNP  -  NUMBER  OF  SYSTEM  NODAL  POINTS 
x  NDOF  -  NUMBER  OF  DEGREES  OF  FREEDOM 

x  N,NN  -  DIMENSIONS  OF  MATRICES  AS  SPECIFIED  IN  DIMENSION 
X  STATEMENTS 

X  E  -  MATERIAL  MODULUS  OF  ELASTICITY 

x  GI  -  SECOND  MOMENT  OF  THE  BEAM  CROSS-SECTION  AREA 

x  ELE  -  ELEMENT  LENGTH 

x  ALPHA  -  LOCATION  OF  ELEMENT  LEFT  GNP 

X  PSIAVE  -  ESTIMATE  OF  PSI 

X  P3ISQ  -  ESTIMATE  OF  PSI  SQUARED 

x  TEND  -  VALUE  OF  TIME  AT  WHICH  THE  SOLUTION  IS  DESIRED 
x  NEQ  -  NO.  OF  FIRST  ORDER  DIFFERENTIAL  EQUATIONS 
x  TIME  -  INDEPENDENT  TIME  VARIABLE 

x  DEL TME  -  TOTAL  TIME  INCREMENT  FOR  ONE  INTEGRATION  STEP 
x  BETA  -  CONSTANT  DETERMINED  BY  BEAM  MATERIAL  PROPERTIES  ONLY 
X  RATE  -  LENGTH  CHANGE  PER  UNIT  TIME 

X  EXEE  -  THE  GLOBAL  NONDIMENSIONAL  AXIS,  THAT  IS,  (X/L) 
x  DELTA  -  THE  VECTOR  OF  NONDIMENSIONAL  NODE  DEFLECTIONS 
X  AND  SLOPES.  MUST  BE  MULTIPLIED  BY  L ( T )  FOR 

x  ACTUAL  DEFLECTIONS.  SLOPES  REMAIN  THE  SAME, 

x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
INCLUDE  'COMMON  FORTRAN' 

DIMENSION  LELTA(NN) , PARAMC NPARAM) ,WKSCNN) 

DIMENSION  YPRIME(NN) 

COMMON  /NORKSP/  RWKSP 

REAL  RWKSPC6608 ) 

EXTERNAL  FCN 
EXTERNAL  FCNJ 
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ODD  O  O  O  O 


0PENC9,  FILE= '/MARK1  INPUT*) 

OPENCIO,  FI L  E=  ' /MARK3B  OUTPUT*) 

0PENC11,  FILE  = '/DATA1  INPUT') 

OPENC 12/  FILE='/DATA2  INPUT') 

CALL  IWKINC6608) 

CALL  ZERO  (DELTA, PARAM) 

CALL  DATA 
CALL  MATX 

PI  =  3 . 1415926  54 

DEFINITION  OF  PARAMETERS  REQUIRED  BY  IMSL  MATH  LIBRARY  ROUTINE 
D IV PAG 

I  DO  =  l 

NEQ=2xNDOF 
TIME-TSTART 
TOL  =  l . OE-4 
PARAMC  4 ) =2000  000 

PARAMC12)  IS  1  FOR  ADAMS  METHOD  AND  2  FOR  GEAR  (STIFF)  METHOD 

PARAM( 1 2  )  =2 
PARAM( 1 3 ) =2 
PARAMC 19 )=1 
PARAM ( 20 ) =NN 

INITIALIZE  THE  DEPENDENT  VARIABLE  ARRAY  DELTA(NEQ)  . 

xx  CAUTION:  THE  NONDIMENSIONAL  VSTAR  IS  CONSTRUCTED  HERE. 

XX  TO  OBTAIN  THE  ACTUAL  INITIAL  DISPLACEMENT  CONFIGURATION,  V, 

XX  SUBSTITUTE  THE  NONDIMENSIONAL  COORDINATE  AXIS  EXEE  IN  THE 
xx  EXPRESSIONS  BELOW  BY  (X/ZLINT)  AND  REPLACE  VSTARxZLINT  BY  V 

IF  (ISTART.EQ.O)  THEN 
IC=NDOF-l 
EXEE  =  0.0 

HR I  TEC  1 0, X)  'THE  INITIAL  TIME  PRIOR  TO  INTEGRATION  =  *,TIME 
HRITEC 10, x)  'THE  INITIAL  DELTA  VECTOR  IS' 

DO  10  1=1, IC, 2 
FI0V2  =  PI/2. 

DELTA(I)  =  0.1  -  0.1XCOSCPIOV2XEXEE) 

DELTACI+1)  =  0 . 1XPI0V2XSINCPI0V2XEXEE) 

DELTA(I)=0.1X(EXEEXX2.0) 

DELTA( 1+1 )  =  0.2XEXEE 

WRITE  (10, X)  DELTA(I) 

WRITE  (10, X)  DELTACI+1) 

EXEE  =  EXEE  +  ELE 
10  CONTINUE 

ELSE 

READ( 1 1 , x)  TSTART.ZL 

WRITE( 1 0 , x )  'RESTART  TIME  =  ' , TSTART,  'WITH  LENGTH  ',ZL 
READC 1 2 , x  )  (DELTACII ) , 1 1  =  1 , NEQ  ) 

WRI TEC  1 0 , x  )  'INITIAL  DELTA  IN  RESTART  FOLLOWS* 

WRITEC 10, x)  (DELTA(JJ), JJ=1,NEQ) 

TIME=TSTART 


END  IF 
IDO  =  l 

CALL  FORMCNEQ, TIME) 


DO  1000  I EQ=1 , NEQ 
YPRIME( IEQ)=0 . 0 
DO  900  IC=1 , NEQ 
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no  o  no 


YPRIME(IEQ)=YPRIME(IEQ)+H(IEQ,IC)XDELTACIC) 

900  CONTINUE 
1000  CONTINUE 

WRITECIO.X)  'INITIAL  YPRIME  VECTOR  PRIOR  TO  ENTRY  TO  DIVPAG' 

WRITE  (10, *)  (YPRIMECIQ),IQ=1,NEQ) 

WRITEC10 ,X)  •  ' 

WRITEC 1 0 , X )  •  • 

WRITECIO.X)  'IN  INTEGRATION  LOOP,  TIME,  LENGTH  AND  DELTA  FOLLOWS' 
WORK  =0.0 

FI  =  RH0*ZLINT*2.*ACC 

VEE  =  6 .x(DELTA(NDOF-3)-DELTA(N DOF-1) )/EL  E2 
&  +  2.*<DELTA(ND0F-2)+2.*DELTA(ND0F))/ELE 
F2  =  (E*GimVEE*x2)/(ZLINT**2) 

F3  =  RHOXC-RATEXDELTA(NDOF)  +  ZL I NTxDEL TA( NEQ- 1 ) ) X*2 
FNEW  =  FI  +  0.5x(F2  -  F3) 

ZL  =  ZLINT 

DO  30  IEND=I , NSTEP 
FOLD  =  FNEW 
ZLOLD  =  ZL 

TEND  =  TSTART  +  DELTMEXFL0AT( IEND  VFLOATC  NSTEP ) 

IFCTEND. GT. 3.0)  GO  TO  35 
CALL  DTIMECIHCUR,  MINUTE,  ISEC) 

I F( I  HOUR  . LT .  18  .AND.  IHOUR  .GE.  7)  GO  TO  35 

CALL  DIVPAG  ( IDO, NEQ, FCN, FCNJ ,G, TIME, TEND, TOL , PAR AM, DELTA) 

IF  ( MODC I  END, 1 ) . EQ . 0 )  THEN 

ZL  =  ZLINT  -  RATEXTIME  +  ACCXCTIMEXX2) 

ZL2  =  ZLX*2 

ZL  DOT  =  -RATE  +  2 . XACCXTIME 
ZLDDOT  =  2.XACC 

ZL  =  9.  +  1.x  COSC  PIxTIME/1.5) 

WRITEC 10, x)  •  ' 

WRITEC 1 0 , X  )  '  ' 

WRITEC 1 0 , x )  'TIME  =  »,  TIME,  'LENGTH  =  ',  ZL 
WRITEC 10,X)  '  ’ 

WRITEC 10,X)  'DELTA  FOLLOWS' 

WRITEC 10, x)  (DELTACIQ), IQ=1,NEQ) 

ZMOM  =  CEXGI/ZL)X(C6./ELE2)X(DELTA(3)-DELTA(1)) 

&  -C2./ELE)X(DELTA(4)+2.XDELTA(2))) 

SHEAR  =  (EXGI/ZL2)X((12./ELE3)X(DELTA(1)-DELTAC3)) 

&  + ( 6 ./ELE2)X(DELTA(2)+DELTA(4))) 

WRITECIO.X)  •  ' 

WRIT  E ( 1 0 , X )  'MOMENT  =  ',  ZMOM,  'SHEAR  =  SHEAR 
END  IF 

FI  =  RH0XZLINTX2 .XACC 

VEE  =  6 . X( DELTAC  ND0F-3 ) -DELTA CNDOF-l ) )/ELE2 
&  +  2.X(DELTA(ND0F-2)+2.XDELTA(ND0F))/ELE 
F2  =  (ExGI)x(VEtx*2)/ZL2 

F  3  =  RHOX(-ZLDOTXDELTACNDOF)  +  ZLXDELT AC NEQ- 1) ) xx2 
FNEW  =  FI  +  0.5XCF2  -  F3) 

DELWORK  =  0 . 5*( FNEW  +  FOLD)X(ZL  -  ZLOLD) 

WORK  =  WORK  +  DELWORK 

WRITECIO.X)  ' ZL  DOT  =  ',  ZLDOT,  'ZLDDOT  =  ', ZLDDOT 
WRITECIO.X)  'OLD  F  =  ' ,  FOLD,  'NEW  F  =  ' ,  FNEW 
WRITECIO.X)  'WORK  =  *,  WORK 

CALL  FORMCNEQ, TIME) 

30  CONTINUE 
35  CONTINUE 

I  D0  =  3 

CALL  DIVPAG  ( I  DO , NEQ , FCN , FCNJ , G , T IME, TEND, TOL , PARAM , DELTA ) 


STOP 

END 


ZERO  ALL  MATRICES 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
SUBROUTINE  ZERO  ( DELTA, PARAMT 

INCLUDE  ‘COMMON  FORTRAN’ 

DIMENSION  DELTA(NN) , PARAM( NPARAM) 

DO  20  1=1, N 
DO  10  J  =  1 ,  N 
AC  I , J ) =0 . 0 
B( I , J  )  =0 . 0 
CCI, J)=0.0 
DC  I , J )  =  0  .  O 
RC I , J )  =  0  .  0 
CONTINUE 
CONTINUE 

DO  AO  1=1, NN 
DEL  TA  C I )  =  0  .  0 
DO  35  J  =  1 , NN 
GC I »  J  )  =0 . 0 
HC I , J  )  =0 , 0 
CONTINUE 
CONTINUE 

DO  50  1=1, NPARAM 
PARAMC I ) =0  .  0 
CONTINUE 

RETURN 

END 

INPUT  DATA 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
SUBROUTINE  DATA 

INCLUDE  'COMMON  FORTRAN' 

READ  C9,X)  NEL , E, GI , RHO, NSTEP, I START , TSTART , ZLINT, DELTME, RATE, ACC 

MR I TEC 6, 'THE  NUMbcR  OF  ELEMENTS  IS  ',  NEL 

NSNP=NEL+1 

ND0F=2XNSNP 


HRITEC 6 , X )  'THE  NUMBER  OF  SYSTEM  NODAL  POINTS  IS  ’,  NSNP 


HRITEC  6 , x ) 
MRITEC  6 , X) 
WRI TEC  6 , x ) 
MRITEC6 , X) 


'THE  NUMBER  OF  DEGREES  OF  FREEDOM  IS 
'THE  MODULUS  OF  ELASTICITY  IS  ',  E 
'THE  MOMENT  OF  INERTIA  IS  ',  GI 
'THE  MASS  PER  UNIT  LENGTH  IS  ',  RHO 


NDOF 


ELE=1 . O/FLOATCNEL ) 

ELE2  =  ELEXX2 
ELE3  =  ELEXX3 
BETA=RHO/CE*GI ) 

MRITEC 6 , x  )  'THE  VALUE  OF  BETA  IS  '  , BETA 


MRITEC  6 , X) 
MRITEC  6 , x  ) 
MRI T  E ( 6 , x ) 
MRITEC  6 , X ) 


'THE  NUMBER  OF  INTEGRATION  STEPS  IS 
'ISTART  IS  1  TCR  RESTART;  HERE  IT  IS 
'THE  INITIAL  LENGTH  IS  ',  ZLINT 
'RATE  OF  AXIAL  MOTION  IS  '  ,  RATE 


NSTEP 

ISTART 


RETURN 

END 


FILL  LARGE  A 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


SUBROUTINE  MATX 


non  o  o 


INCLUDE  'COMMON  FORTRAN' 

C  CALCULATE  LITTLE  A  MATRIX 

A11=12.0/(ELE**3.0) 

A1 2=6 .0/CELEXX2.0) 

A1 3= ( -1 . C  )XA1 1 

A14=A12 

A21 =A12 

A22  =  4 . 0/ELE 

A23= ( -1 . 0  )*A12 

A24  =  2 . 0/ELE 

A31=A13 

A32=A23 

A33=A 1 1 

A  3  4  =  A  2  3 

A41-A14 

A42=A24 

A43=A34 

A44  =A22 

C  CALCULATE  THE  ELEMENTAL  B  MATRIX 

B 1 1  =  ( -6 . 0)/(5. OXELE) 

BI2  =  -1  .1 
B13  =  ( -1  .  0)XB11 
B1 4  =- . 1 
B2 1  =  B 1 4 

B22=(-2. 0*ELE)/15 . 0 
B23=  .  1 

B24=ELE/30.0 
B  31  =  B1 3 
B32=B23 
B33=BI 1 
B34  =  l  .  1 
B41 =B21 
B42=B24 
B43= . 1 
B44=B22 

CALCULATE  THE  ELEMENTAL  C  MATRIX 


Cll  =  -  .  5 
C12=ELE/10.0 
Cl  3= . 5 

C14=(-.1)XELE 
C21 =C1 4 
C22  =  0 . 0 
C23=C12 

C24  =  ( ELE**2)/(-60 . 0  ) 

C31  =  C1 1 
C32=C21 
C33=C1 3 
C34=C1 2 
C  4  I  =  C 1  2 

C42=C-1 . 0)XC24 
C  4  3  =  C  3  2 
C44 =C22 

CALCULATE  THE  ELEMENTAL  D  MATRIX 

Dll=13.0XELE/35.0 
D12=(11,0/210.0)X(ELE*X2.0) 
D13  =  9 . 0XELE/70  .  0 
Cl  4  =  C -13. 0/420 . 0)*( CL  Ex *2) 
D2 1 =  D1 2 

D22  =  (ELEXX5.0)/(105.  ) 

D23=(13.0XELExx2)/420. 

D24=-(ELEXX3.0)/(140.0) 

D  3 1  =  D 1 3 
D32=D23 
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D  3  3  =  D 1 1 

D34  =  (-11.0*ELE**2)/'210. 

D41 =D1 4 
D42=D24 
D43=D34 
D44=D22 

C  FILL  THE  GLOBAL  A,B,C,  AND  D  MATRICES 
LNEL  =  2#NEL -1 
DO  10  1=1,  LNEL,  2 
ALPHA  =  0  .  0 

PSIAVE=ALPHA+(ELE/2.Q) 

P3ISQ=PSIAVE**2 

ACI,I)=ACI,I)+A11 
A(I,I+1)=A(I,I+1)+A12 
AC  I , I+2)=A13 
AC  1 , 1+5 ) =A1 4 

AC  1  +  1 , I  )=A( 1  +  1 , I  )  +  A21 
AC  1  +  1 , 1  +  1  )=A( 1  +  1 , 1  +  1 )+A22 
AC I+I , I+2)=A23 
AC  1  + 1 , 1  +  3 ) =A24 

AC  1+2, 1 )=A31 
AC  1  +  2, 1  +  1 )  =  A32 
AC  1  +  2, 1  +  2  )rA33 
A(I+2,I+3)=A34 

AC  1  +  3 , 1 ) =A41 
AC  1  +  3, I  +  I  )=A42 
AC  1  +  3, I  +  2)=A43 
AC  1  +  3, I  +  3)=A44 


BCI,I)=B(I,I)+(B11*P$ISQ) 

BC 1 , 1  +  1 ) =BC I , I  +  I  )  +  C  B12*PSISQ  ) 

B(I,I+2)=B13*P$ISQ 

BCI,I+3)=B14*PSISQ 

BCI+1,I)=B(I+1,I)+(B21*PSISQ) 

BC 1  +  1, 1  +  1  )  =  B( 1  +  1, 1  +  1  )  +  CB22*PSISQ) 

B(I+1,I+2)=B23*PSISQ 

BCI+1, I+3)=B24*PSISQ 

BC 1  +  2, 1  )  =  B31#PSISQ 
BC 1+2, 1+1 ) =  B32*PS I SQ 
BC 1+2, I+2)=B33*PSIS9 
BC 1+2, I+3)=B34*PSISQ 

BCI+3,I)=B41*PSISQ 
BC  1  +  3, 1  +  1 )=B42*PSISQ 
BC 1+3, I+2)=B43*PSISQ 
BCI+3,I+3)=B44*PSI5Q 

C(I, I)=C(I , I )+(Cll*PSIAVE) 

CCI ,  I  +  I  )=CC  1 , 1  +  1  > +  (  Cl  2*P3IAVE) 

CC I , I+2)=C13*PSIAVE 
CCI,I+3)=C14*PSIAVE 

CC 1+1 , I ) =CC I+I , I)+CC21*PSIAVE) 
CCI+1, I+1)=C(I+1, 1+1 )+(C22*P3IAVE) 
CCI+1, 1+2 ) =C25*PSIAVE 
CCI+1, I+3)=C24*PSIAVE 

CC 1  +  2, 1  )=C31*PSIAVE 
CC 1  +  2, 1  +  1  )=C32*P3IAVE 
CCI+2, I+2)=C33*PSIAVE 
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nnn 


C( 1+2, I+3)=C34*PSIAVE 


CCI+3, I)=C41*PSIAVE 
C  C  1  +  3, I  +  1 )=C42*PSIAVE 
CCI+5,I+2)=C43*PSIAVE 
C(  1  +  3, 1+3 ) =C44*PSIAVE 

DCI,I)=D(I,I)+D11 

D(I,I+l)=D(I,I+lJ+D12 

DCI,I+2)=D13 

D(I,I+3)=D14 

D(I+1,I)=DCI+1,I)+D21 
D( 1+1 , 1+1 )  =  D( 1+1 , 1+1 )+D22 
DC  1  + 1 , 1  +  2 ) =D23 
DC  1  +  1 , 1+3 )-D24 

DC  1  +  2, 1 )=D31 
DC  1+2, 1  +  1 ) =  D32 
DC  1+2, 1  +  2 1  =  D  3  3 
DC  1  +  2, 1  +  3 )  =  D34 

DC  1  +  3, 1 )  =  D41 
DC  1  +  3, 1  +  1 )  =  D42 
DCI+3,I+2)=D43 
DC  1  +  3, I+3)=D44 

ALPHA=ALPHA+ELE 

10  CONTINUE 

RETURN 

END 


x*xxxxx¥xxxx»¥x*xxxxxx*#xxx*xxxxxxxxxxxxxxxxxx*xx#xxx*x#***xxxxx 
SUBROUTINE  FORMC  NEQ , TIME ) 

INCLUDE  'COMMON  FORTRAN' 

ZL  =  ZLINT  -  RATE*TIME  +  ACC*( 7IME**2 ) 

ZL  DOT  =  -RATE  +  2.*ACC*TIME 
ZLDDOT  =  2 . *ACC 
FI  =  3.141592654 
PARAM  =  PI/1 .5 

ZL  =  9 .  +  1 . *  COSC  PARAM#TIME ) 

Z L  DOT  =  - C  PARAM ) *SI NC  PARAMxTIME ) 

ZLDDOT  =  -CPARAM**2)*C0SC  P ARAM* T I ME) 

AC0EFF=C-1 . 0)/( CZL«*4)*BETA) 

ECOEF F -  -( CZLDOT/ZL )«*2) 

CCGEFF=  -C2 . *CZLD0T/ZL  )**2)  +  ZLDDOT 

DO  20  I = 1 , NDOF 
CO  15  J=l,NDOF 

R ( I . J  )  -  A  (  I ,  J  ) *ACOEFF  +  B C I , J  )  * BCOEFF  +  CC I , J ) XCCOEFF 
CCCI, J)=  -2.*BETA*CZLD0T/ZL)*CCI, J) 

15  CONTINUE 

20  CONTINUE 

C  REDUCE  SYSTEM  OF  EQUATIONS  FROM  SECOND  TO  FIRST  O.D.E. 

DO  100  1=1, NDOF 
GC I , I  )  =  1  .  0 
100  CONTINUE 

DO  300  1=1, NDOF 
K=I+NDOF 
DO  200  J  =  1 , NDOF 
L  L  = J  +  NDOF 
GC  K , L  L )  =  D( I , J ) 

2  C  0  CONTINUE 
300  CONTINUE 
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DO  <400  1  =  1  ,  NDOF 
1 1  =  I +  NDOF 
H! 1 , 1 1 )  =1 . 0 
400  CONTINUE 

DO  600  I=i , NDOF 
K= I +NDOF 

DO  500  J  =  1 , NDOF 
HIK, J)=R! I, J) 

500  CONTINUE 

600  CONTINUE 

DO  700  1=1, NDOF 
I I = I+NDOF 

DO  650  J  =  1 , DOF 
J J  =  J  +  NDOF 

HCII.JJ)  =  CCCI.J) 

650  CONTINUE 
700  CONTINUE 

C  IMPOSE  FIXED  END  BOUNDARY  CONDITIONS 

DO  800  J  =  1 , 2 

DO  750  1=1, NEQ 

G( J  +  NDOF , I )  =  0 . 0 
H( J , I  )  =0 . 0 
H ! J  +  NDOF , I ) =0 . 0 
750  CONTINUE 
800  CONTINUE 

G( NDOF+1 , NDOF+1 ) =1 . 0 
G!ND0F+2,NDQF+2)=1.0 

RETURN 

END 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

SUBROUTINE  FCN  < NEQ, TIME, DELTA, YPRIME) 

INCLUDE  'COMMON  FORTRAN’ 

DIMENSION  YPRIME! NEQ), DELTA C NEQ) 

REAL  L 

DELTA!  1  )  =  0.0 
DEL  T  A ( 2  )  =  0.0 
DEL  T  A( NDOF+1 )  =  0.0 
DELTA<ND0F+2)  =  0.0 

C  FORM  YPRIME 

DO  1000  IEQ=1 , NEQ 
YPRIME! IEQ)=0 . 0 
DO  900  I C  =  1 , NEQ 

YPRIME!IEQ)=YPRIME!IEQ)+H!IE<3,IC)*DELTA!IC) 

900  CONTINUE 
1000  CONTINUE 

RETURN 

END 

FUNCTION  FCNJINEQ,  TIME,  DELTA,  PD) 

REAL  'IME,  DELTA(NEQ),  PD(x) 

fcnj  =  o  .  o 

RETURN 

END 
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